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The  work  carried  out  \mder  this  contract  may  be  subdivided  according  to  the  following 
topics: 

1)  Grid  Generation  and  Pre-Processing; 

2)  Flow  Solvers; 

3)  Grid  Adaptation; 

4)  Rigid  Body  Motion; 

5)  Visualization  and  Animation; 

6)  Demonstration  Calculation  and  Results. 

The  description  of  the  main  accomplishments  of  the  present  work  are  listed  according 
to  these  topics  in  the  following. 

1.  GRID  GENERATION  AND  PRE-PROCESSING 

All  the  unstructured  grids  required  for  the  simulation  of  compressible  flows  with  moving 
bodies  were  generated  using  the  advancing  front  technique  [1].  At  the  beginning  of 
the  present  contract,  we  bad  a  working,  but  slow  and  unreliable  unstnictxired  grid 
generator.  The  fifllowing  improvements  to  this  original  capability  were  implemented  as 
part  of  this  grant. 

1.1  Reliable  3-D  Grid  Unstructured  Generator:  We  improved  the  reliability  of  the  3-D 
grid  generate  to  a  point  where  it  could  be  used  for  the  many  regenerations  that  typi¬ 
cally  take  place  during  a  transient  flow  simulation  with  moving  bodies.  This  significant 
increase  in  reliability  was  achieved  by: 

a)  not  allowing  any  bad  elements  during  the  generation  process;  and 

b)  enlarging  and  remeshing  those  regions  where  new  elements  could  not  be  introduced. 
Thus,  we  first  attempt  to  complete  the  mesh,  skipping  those  faces  that  do  not  give  rise 
to  good  elements.  If  pockets  of  unmeshed  regions  remain,  we  enlarge  them  somewhat, 
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and  regrid  them.  This  ‘sweep  and  retry’  technique  has  proven  extremely  robust  and 
reliable.  In  fact,  the  3-D  unstructured  grid  generator  has  not  failed  a  sin^e  time  since 
this  capability  was  implemented. 

1-2  Incorporation  of  Stretching  in  3-D  Grid  Generator:  The  ability  to  generate 
stretched  elements  has  the  potential  payoff  of  order-of-magnitude  reductions  in  the 
required  number  of  elements  and  CPU  time.  Stretching  wais  taken  into  account  during 
all  phases  of  the  grid  generation  process:  when  generating  sides  along  lines,  triangles 
on  the  siirfaces,  and  tetrahedra  in  the  domain.  The  most  difficult  of  these  was  the 
generation  of  triangles  on  surfaces. 

1.3  Fast  3-D  Grid  Generator:  A  typical  CFD  run  with  moving  bodies  requires  many 
regriddings.  Thus,  it  may  cost  as  much  to  advance  the  flow  solver  as  to  regrid.  There¬ 
fore,  the  grid  generator  has  to  be  made  as  fast  as  possible.  The  following  techniques 
were  used  to  improve  the  performance  of  the  advancing  front  grid  generator: 

a)  Filtering:  Tjrpically,  the  number  of  close  points  and  faces  is  fw  too  conservative, 
i.e.  large.  As  an  example,  consider  the  search  for  close  points:  there  may  be  up  to 
eight  points  inside  an  octant,  but  of  these  only  one  may  be  close  to  the  face  to  be 
taken  out.  The  Mea  is  to  filter  out  these  ‘distant’  faces  and  points  in  order  to  avoid 
extra  work  afterwards.  While  the  search  operations  are  diffinilt  to  vectorize,  these 
filtering  operations  lend  themselves  to  vectorization  in  a  straightforward  way,  leading 
to  a  considerable  overall  reduction  in  CPU  requirements. 

b)  Automatic  Reduction  of  Unused  Points:  As  the  front  advances  into  the  domain  and 
more  and  more  tetrahedra  are  generated,  the  number  of  tree-levels  increases.  This 
automatically  implies  an  increase  in  CPU-time,  as  more  steps  are  required  to  reach  the 
lower  levels  of  the  trees.  In  order  to  reduce  this  CPU-increase  as  much  as  possible, 
all  trees  are  automatically  restructured.  All  points  which  eire  completely  surrounded 
by  tetrahedra  are  eliminated  from  the  trees.  This  procedure  has  proven  to  be  ex¬ 
tremely  effective.  It  reduces  the  asymptotic  complexity  of  the  grid  generator  to  less 
than  O(JV’logJV).  In  fact,  in  most  practical  cases  one  observes  a  linesu:  0{N)  asymp¬ 
totic  complexity,  as  CPU  is  traded  between  subroutine  call  overheads  etnd  less  close 
faces  on  average  for  large  problems. 

c)  Global  H-refinement:  While  the  basic  2uivancing  front  algorithm  is  a  scsdar  algo¬ 
rithm,  h-refinement  can  be  completely  vectorized.  Therefore,  the  grid  generation  pro¬ 
cess  can  be  made  considerably  faster  by  first  generating  a  coarser,  but  stretched  mesh, 
smd  then  refining  globally  this  first  mesh  with  classic  h-refinement  [2].  Typical  speed- 
ups  achieved  by  using  this  approach  are  1:6  to  1:7. 

1.4  Local  Remeshiny:  Practical  simulations  reveeded  that  the  appearance  of  badly  dis¬ 
torted  elements  occurred  at  a  frequency  that  was  much  higher  than  expected  from  the 
element  size  prescribed.  Given  the  relatively  high  cost  of  global  remeshing,  we  explored 
the  idea  of  local  remeshing  in  the  vicinity  of  the  elements  that  became  too  distorted. 
Thus,  wherever  the  elements  become  too  distorted,  we  open  up  ‘holes’  in  the  mesh.  We 
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then  recompute  the  error  indicators,  and  remesh  adaptively  the  ‘holes’  using  the  ad¬ 
vancing  front  method.  Typically,  only  a  very  small  number  of  elements  (<  10)  becomes 
so  distorted  that  a  remeshing  is  required.  Thus,  local  remeshing  is  a  very  economical 
tool  that  has  allowed  us  to  reduce  CPU-requirements  by  more  than  60%  for  typical 
runs. 

1.5  Mesh  Smoothing  and  Ontimizatinn;  After  the  advancing  grid  generator  has  filled 
the  region  to  be  meshed  with  elements,  it  is  advisable  to  smooth  the  mesh  in  order  to 
improve  the  element  quality  further.  This  not  only  avoids  bad  results  due  to  deformed 
elements,  but  also  allows  larger  timesteps  to  be  taken,  significantly  reducing  CPU 
requirements.  Two  different  ways  of  optimizing  or  smoothing  unstructured  grids  were 
explored. 

1.5.1  Spring  System  Analogy  with  Local  Regeneration:  The  spring  system  analogy  has 
been  widely  used  for  smoothing  unstructured  grids.  Each  side  or  edge  in  the  mesh  is 
supposed  to  represent  a  spring.  Thus,  the  force  acting  on  each  point  is  given  by: 

fj=C^(x,  -X.)  ,  (1) 

J=1 

where  c  denotes  the  spring  constant,  Xi  the  coordinates  of  the  point,  and  the  sum 
extends  over  all  siurotmding  points.  The  time-advancement  for  the  coordinates  is 
accomplished  as  follows: 

Axj  =  At~f,  .  (2) 

ns, 

At  the  boundary  of  the  domain.  Ax  =  0.  Usually,  3-4  passes  over  the  mesh  yield  an 
acceptable  mesh.  Unfortunately,  this  type  of  smoothing  can  produce  negative  elements. 
Typically,  these  axe  very  few  compared  to  the  overall  number  of  elements  (<  1%).  This 
makes  local  remeshing  an  attractive  option.  Therefore,  we  treat  these  negative  elements 
as  distorted,  remove  them,  and  regrid  these  'boles’  again.  To  our  knowledge,  no  other 
grid  generator  presently  in  operation  uses  smoothing  on  a  routine  basis.  This  has 
enabled  us  to  produce  the  best  imstructured  grids  presently  available,  as  the  statistics 
clearly  show.  Figure  1  is  an  example  of  such  a  statistic. 

1.5.2  Smoothing  vi^  The  idea  here  is  to  state  a  functional  that  describes 

the  optimal  mesh.  Then,  the  present  mesh  is  optimized  using  a  PoUak-Ribiere  conjugate 
gradient  algorithm.  This  technique  has  been  used  extensively,  and  very  successfully,  for 
the  optimization  of  structured  grids.  However,  unlike  most  finite  difference  techniques, 
it  is  completely  general,  and  thus  can  be  extended  to  unstructured  grids.  Its  main 
advantage  is  that  it  will  never  give  rise  to  elements  with  negative  Jacobians.  This 
technique  has  been  described  in  detail  in  Refs. [3-5].  Starting  from  a  mesh  that  was 
already  smoothed  using  the  spring  system  analogy,  one  observes  always  an  increase 
in  quality.  In  some  cases,  this  improvement  is  dramatic,  leading  to  increases  in  the 
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Fill  .  Mall  Quality  Statistic.  This  Figure  Demonstrates  the  Significant  Mesh-Quaiity 
Improvement  Obtained  using  the  Advanced  3-D  Smoothing  Algorithm  Recently  Developed 


allowable  timestep  of  an  order  of  magnitude  [5].  We  have  included  more  information 
on  this  optimization  technique  in  Appendix  1-4. 

1.6  Mesh  Post-Processor:  Practical  runs  show  that  even  after  smoothing  or  otherwise 
improving  the  mesh  a  few  badly  distorted  elements  may  still  remain.  Remember  that 
in  3-D,  where  even  grids  for  invisdd  flow  calculations  exceed  1  million  tetrahedra,  every 
small  possibility  becomes  a  reality.  These  badly  deformed  elements  that  remsdn  at  the 
end  are  removed  with  a  special  post-processor  [5].  The  elements  in  question  £ire  deleted 
by  collapsing  one  of  their  sides  and  removing  one  of  the  points  without  creating  any 
elements  of  worse  shape. 

1.7  Diagnostics:  An  extensive  smte  of  diagnostics  subroutines  that  check  for  input 
errors  was  added  to  the  grid  generator.  The  most  common  types  of  errors  (double 
point  definition,  unclosed  surfaces,  surface  orientation,  stirface  crossings,  incomplete 
boundary  condition  definition,  etc.)  are  checked  internally  by  the  grid  generator  at  nm- 
time.  Our  expferience  has  been  that  these  diagnostics  not  only  significsmtly  improves 
the  reliability  of  the  CFD  nms,  but  also  reduces  problem  set-up  times  considerably  by 
pointing  the  user  to  potential  problems. 

1.8  PREGEN /PREBACK:  During  the  course  of  the  present  effort  it  became  clear 
that  the  error-free  problem  setup  was  the  most  real-time  intensive  factor  of  a  typical 
CFD  run  with  moving  bodies.  Therefore,  a  considerable  amount  of  effort  was  spent  in 
reducing  the  burden  of  specifying  the  problem  to  be  solved.  Two  tools  were  developed 
to  this  end: 

1.8.1  PREGEN:  PREGEN  consists  of  a  suite  of  subroutines  that  allows  the  user  to 
produce  grid  generator  compatible,  error-free  input  in  a  faster  way.  PREGEN  not 
only  allows  the  user  to  exercise  basic  CAD-CAM  operations  (shrinking,  translations, 
rotations,  surface  lofting,  etc.),  but  also  eases  the  merging  of  several  parts  of  the  surface 
into  one  cohesive,  well-defined  input-file.  This  allows  the  merger  of  files  produced  by 
different  users  and/or  different  stirface  generators.  PREGEN  has  a  whole  series  of  built- 
in  diagnostics  to  avoid  such  undesirable  features  as  doubly  defined  points,  isolated 
points  or  lines,  badly  defined  lines  or  surfaces,  etc.  PREGEN  is  highly  interactive, 
allowing  the  user  to  monitor  every  step  during  problem  set-up  in  a  graphical  way. 

1.8.2  PREBACK:  PREBACK  generates  3-D  background  grids  by  lofting  or  rotating 
a  2-D  triangulation.  In  addition,  it  allows  interactive  operations  to  move  background 
grid-points  around  (translation,  rotation,  shrinking,  etc.),  as  well  as  to  modify  the  grid- 
generation  parameters  in  space  (size,  shape).  By  being  able  to  see  both  the  surface¬ 
defining  CAD-CAM  data  and  the  backgrotmd  grid,  the  user  can  quickly  generate  a 
background  grid  that  is  suited  for  the  desired  distribution  of  element  size,  stretching 
and  stretching  directions  in  3-D  space. 


2.  FLOW  SOLVERS 


All  the  flow  solvers  used  for  the  present  effort  are  based  on  finite  elements  and  operate 
on  unstructured  grids.  At  the  beginning  of  the  period  of  performance,  the  flow  solvers 
were  based  on  FEM-FCT  techniques  [6],  and  worked  on  element  loops.  This  type  of 
flow  code  has  been  the  workhorse  for  many  of  the  runs  done.  In  the  last  hzdf  year 
we  have  developed  some  new  flow  solvers  that  were  aimed  specifically  at  reducing  the 
indirect  addressing  (i/a)  cost  and  the  memory  reqmrements.  These  newer  flow  codes 
are  based  on  an  edge  data-structure,  and  use  new  renumbering  techniques  to  reduce 
both  the  cache-misses  as  well  as  to  increase  the  performance  on  vector  machines.  The 
main  accomplishments  in  this  area  are  detailed  in  the  following. 

2.1  ALE  Capability:  In  order  to  handle  the  moving  frames  of  reference  associated  with 
the  moving  finite  elements,  the  partial  differential  equations  need  to  be  modified.  This 
is  most  easily  accomplished  by  the  Arbitrary  Lagrangian-Eulerian  (ALE)  formulation. 
The  derivation  of  the  equations  may  be  foimd  in  [7].  Here,  we  just  state  the  final  form 
of  the  equations  of  motion.  Given  the  velocity  field  w  for  the  elements 

w  =  (ttf*,u;*,u>*)  ,  (3) 

the  Euler  equations  that  describe  an  inviscid,  compressible  fluid  may  be  written  as 
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Observe  that  in  the  case  of  no  element  movement  ( w  =  0),  we  recover  the  usual  Eulerizin 
conservation-law  form  of  the  Euler  equations.  If,  however,  the  elements  move  with  the 
particle  velocity  (w  =  v),  we  recover  the  Lagrangian  form  of  the  equations  of  motion. 
From  the  numerical  p<^t  of  view,  E)qn.(4)  implies  that  sdl  that  is  required  when  going 
from  an  Euleriaa  frame  to  an  ALE-frame  is  a  modified  evaluation  of  the  fluxes  on  the 
left-hand  side,  and  the  additional  evaluation  of  soxirce- terms  on  the  right-hand  side. 
These  ALE  equations  were  implemented  in  all  the  flow  solvers  used  for  the  present 
effort. 

2.2  Edge-Based  Solvers:  A  significant  reduction  in  indirect  addressing  (i/a)  costs  can 
be  realized  by  going  from  an  element  to  an  edge-based  data  structure  for  the  flow 
solver.  We  developed  a  new  series  of  edge-based  flow  solvers.  When  developing  these 
new  codes,  we  incorporated  all  the  coding  lessons  that  we  learned  over  the  years.  The 
latest  code,  FEFL093,  runs  at  a  sustained  rate  of  115  MFlops  on  the  CRAY-YMP,  and 
has  significantly  less  Flops  per  update  than  the  previous  code  (FEFL054).  At  the  same 
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time,  i/a  costs  were  reduced  by  a  factor  of  3.5.  We  plan  to  continue  the  development 
of  this  new  c^ability  in  the  future,  and  expect  to  reduce  i/a  costs  by  another  factor 
of  3  shortly. 

2.3  Renumbering  Stratesnes;  This  development  was  carried  out  specifically  for  machines 
with  small  cache-memory,  like  the  CONVEX  or  the  INTEL  hypercubes.  The  idea  is  to 
still  be  able  to  operate  in  vector-mode,  but  to  avoid  large  jumps  in  data  accessed.  Of 
the  many  techniques  explored,  the  following  two  were  the  mosts  successful: 

2.3.1  Element  Renumbering  According  to  Smallest  Point:  The  elements  are  renum¬ 
bered  according  to  their  smallest  point.  This  technique  by  itself  reduced  CPU  times 
for  typical  runs  by  50%  on  the  IBM  RISC-6000/530  for  large  grids. 

2.3.2  Small  Group  Element  Colouring:  In  order  to  avoid  memory  conflicts  inside  vector- 
loops,  elements  or  edges  have  to  be  grouped  together,  or  ‘coloured’.  For  a  CRAY,  the 
aim  is  to  achieve  as  small  a  number  of  groups  with  as  large  a  vector  as  possible.  For 
a  CONVEX,  such  a  modus  operand!  will  lead  to  a  large  number  of  cache-misses.  The 
idea  is  to  reduce  the  vector  lengths  to  64  or  128,  but  to  operate  on  local  data  as  much 
as  possible.  This  is  accomplished  by  first  renumbering  according  to  smallest  point,  eind 
then  by  colouring  in  small  groups,  always  starting  at  the  beginning. 

2.4  Linelets:  Any  implicit  flow  solver  (and  we  will  need  these  for  Navier-Stokes  simu¬ 
lations)  requires  the  solution  of  a  large  system  of  linear  equations  of  the  form 

Ku  =  r  ,  (5) 

where  K  is  the  matrix  arising  from  implicit  timestepping,  u  the  desired  vector  of 
unknowns,  and  r  the  right  hand  side  vector.  The  fzatest  way  to  solve  such  large  sys¬ 
tems  of  equations  is  through  the  use  of  unstructured  multigrid  solvers.  They  require 
good  smoothers,  as  well  as  efficient  intergrid  transfer  operators.  During  the  present 
year,  we  developed  a  class  of  iterative  solvers  that  lie  between  the  complexity  of  mul¬ 
tiple  grids  and  the  excessive  memory  reqxiirements  of  direct  solvers.  Codes  based  on 
structured  grids  have  explored  for  a  long  time  the  useful  properties  of  line-relaxation. 
Line-relaxation  offers  an  economical  way  to  circumvent  directional  stifihess  by  joining 
together  neighbors  of  neighbors  along  the  line.  Thus,  it  offers  a  practical  way  to  derive 
good  preconditiooers  and  smoothers.  The  concept  of  lines  translates  to  snakes  in  the 
context  of  an  unstructured  grid.  The  corresponding  relaxation  scheme  to  solve  Eqn.(5) 
becomes 


KrIC2Au  =  r  —  Ku  .  (6) 

Here  Ki,  Kj  denote  the  entries  of  K  for  the  active  point-point  combinations  that  define 
the  snake.  Observe  that  this  formulation  is  not  dimensionally  consistent.  Therefore,  a 
better  formulation  is  to  use 


K  =  Di(I  -h  E)D*  «  0^(1  +  E,  )(I  -H  E2)Di  , 


(7) 


Figure  1:  Snake  and  Linelets  for  a  2-D  Discretization 
Preferred  Direction  ( i ,0) 


where  D,E  denote  the  diagonal  and  ofF-diagonal  entries  of  K.  Thus,  the  relaxation 
scheme  becomes 

D*(I  +  Ei)(I  +  E2)D^Au  =  r-Ku  .  (8) 

In  order  to  achieve  higher  rates  of  convergence,  we  use  Eqn.(8)  within  a  preconditioned 
Conjugate  gradient  algorithm.  As  an  xxnstructured  grid  usually  does  not  possess  an 
equal  ntunber  of  gridpoints  along  a  certain  direction,  the  resulting  snakes  may  often 
exhibit  folding  (see  Figure  2b).  This  implies  that  the  information  flow  from  the  do¬ 
main  to  the  botmdary  may  be  slowed  down  considerably.  This  is  not  important  for 
smoothers,  but  crucial  for  the  preconditioners  required  in  one-grid  solvers.  In  order  to 
obtain  a  steady  flow  of  information  towards  the  boundaries  in  all  directions,  we  recon¬ 
nect  the  snake  in  the  direction  it  intended  to  continue  wherever  it  folds.  This  gives 
rise  to  a  more  complex  structure,  which  we  call  linelet  (see  Figiire  2c).  Whereas  the 
storage  requirements  of  snakes  are  fixed  (3iV,  where  N  is  the  niunber  of  unknowns),  the 
storage  requirements  of  linelets  depend  on  the  structure  of  the  mesh  and  the  renum¬ 
bering  chosen.  Using  reverse  Cuthill-McKee  ordering,  as  well  as  octrees,  one  observes 
0(5iV  —  lOiV)  storage  requirements.  This  is  deemed  2u:ceptable,  as  it  takes  much  more 
than  6N  operations  to  build  a  new  right-hand  side.  As  expected,  the  convergence  rate 
of  the  preconditioned  Conjugate  gradient  algorithm  increases  considerably  when  going 
from  snakes  to  linelets.  This  completely  new  concept  went  through  several  stages  of 
development.  Starting  with  snakes,  we  soon  realized  the  problems  with  fol<^g,  amd 
had  to  develop  ways  to  construct  the  linelets  in  an  efiicient  way.  Because  the  size  of 
the  system  of  equations  to  be  solved  depends  on  the  numbering  of  nodes,  near-optimal 
renumbering  schemes  had  to  be  tested.  We  went  through  four  generations  of  renumber¬ 
ing  strategies.  Compared  to  our  current  scheme,  a  simple-minded  renumbering  would 
yield  equation  systems  that  require  4-8  times  more  storage  and  CPU  for  the  same  per¬ 
formance.  So  far,  we  have  used  this  linelet-based  solver  for  scalar  elliptic  problems.  In 
the  coming  year,  we  will  extend  it  to  the  fully  coupled  compressible  Navier-Stokes  case. 

2.4.1  Vectorization  of  Linelet-Matrix  Solution:  The  first  linelet-matrix  solvers  we  im¬ 
plemented  were  based  on  the  general  Crout  LU  decomposition,  which  can  be  found 
in  many  Finite  Element  textbooks.  Performance  tracing  on  the  CRAY-YMP  showed 
that  this  portion  of  the  code  was  running  at  about  2.5Mflops,  dismal  for  a  machine 
that  is  rated  at  2S0Mflop8  per  processor.  This  veiy  low  speed  was  the  result  of  short 
vector  lengths  in  the  mainly  tridiagonal  structure  of  the  linelet-matrix.  Vectorization 
of  both  the  LU  factorization,  as  well  as  the  forward  and  backward  reduction  during 
solution  can  be  accomplished  if  one  processes  in  parallel  as  many  linelets  as  possible. 
The  entries  in  the  linelet-matrix  are  categorized  into  diagonal  (points  with  prescribed 
unknowns),  tridiagonal  or  non-tridiagonal.  At  each  stage,  all  active  non-tridiagoned 
inhibitors  are  processed  first.  Thereafter,  all  available  tridiagonal  entries  are  processed 
in  vector-mode.  Towards  the  end  of  the  factorization  or  reduction  process,  the  available 
number  of  tridiagonal  locations  may  diminish  to  a  point  where  scalar  processing  of  the 
remaining  entries  is  faster  than  extremely  short  vector- loops.  For  this  reason,  a  scalar 
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tridiagonal  category  is  added  to  the  complete  solver.  The  complete  factorization  or 
reduction  process,  illtistrated  in  Figure  3,  is  then  given  by: 


1.  Identify  ‘Origins’  (V) 

2.  Solve  For  Non-Tridiag.  ‘Inhibitors’  (SV,DOT) 

3.  Solve  For  Available  Tridiagonals  (V) 

4.  Solve  For  Remaining  Linelets  (FS) 

5.  If  Not  Finished:  GOTO  2. 


The  vectorized  version  of  the  linelet-solver  presently  runs  at  ISMflops.  This  represents 
a  speed-up  of  7  as  compared  to  the  scalar  version,  but  is  still  deemed  unacceptable. 
FWther  research  is  being  devoted  to  this  subject  at  the  present  time.  Of  covuse,  for 
the  solution  of  blocks  of  unknowns  (as  required  for  the  compressible  Navier-Stokes 
equations),  as  compared  to  a  single  variable,  the  achievable  Mflop-rate  is  expected  to 
be  significantly  higher. 

2.5  Coding  of  Body  Force  Reduction  Subroutines:  In  order  to  couple  the  flowfields  with 
the  rigid  body  motion  one  has  to  extract  the  forces  and  moments  exerted  by  the  fluid 
on  the  bodies  present.  These  data  reduction  subroutines  were  coded  and  tested. 

2.6  Turbulence  Model:  The  algebraic  Baldwin-Lomax  turbulence  model  [8]  was  imple¬ 
mented  in  the  2-D  flow  solver.  This  is  the  simplest  way  to  add  the  effects  of  turbvilence. 
The  model  is  straightforward  to  implement  as  the  only  difference  in  the  flow  equations  is 
the  increase  of  viscosity.  The  model  represents  the  inner  and  outer  parts  of  a  boimdary 
layer  as  follows: 

Inner  part: 

Ht  p\ui\  [fcy(l  -  exp{-y/A))]^ 

Outer  part: 

Pt  ~  rnaxl  > 

with 


F(y)  =  y|u;|(l-exp(-y/A))  . 

The  exponential  factor  in  the  inner  part  is  the  V2ui  Driest  damping  factor  which  matches 
the  damping  of  the  wall.  In  the  second  equation,  Fmat  is  computed  along  a  normal 
to  the  wall.  This  is  particularly  easy  when  dealing  with  structured  grids,  but  requires 
some  effort  for  xinstructured  grids.  The  required  data  structures  were  discussed  by 
Rostand  [9].  One  complete  evaluation  of  the  turbulent  viscosity  requires  the  following 
transfer  of  information: 

a)  Vorticity  from  elements  or  points  in  the  mesh  to  the  appropriate  normals  to  the  walls. 
Along  each  normal,  the  vorticity  |u;|  is  required  to  evaluate  the  turbulent  viscosity.  This 


tr2uisfer  of  vortidty  is  accomplished  with  a  linked  list  of  intersections  of  normals  to  the 
wall  with  elements  of  the  mesh.  This  list  is  constructed  by  starting  from  the  surface  and 
moving  along  the  normal.  All  that  is  required  to  do  so  is  a  Ust  of  elements  surrovinding 
elements. 

b)  Transfer  of  the  turbulent  viscosity  from  the  normals  to  the  walls  to  the  elements  or 
points  of  the  mesh.  In  the  present  case,  we  transfer  to  points.  This  not  only  reduces  the 
transfers  required  (there  are  less  points  than  elements  in  a  mesh),  but  also  introduces  a 
beneficial  smoothing  effect  at  element  level.  For  each  point  in  the  mesh  close  to  wetted 
surfaM:es,  we  find  the  closest  normals  surrounding  it,  and  then  the  two  closest  points 
along  each  of  the  normals.  Thus,  a  point  in  the  mesh  assembles  the  turbulent  viscosity 
from  four  points  along  the  normads  of  a  wetted  surface.  The  search  for  the  closest  points 
is  done  using  quad-trees  [10].  In  the  case  of  the  mixing  of  several  boundary  layers,  /if 
is  computed  from  the  contribution  of  each  wall  weighted  by  its  distance  to  the  point. 

3.  ADAPTIVE  REFINEMENT 

Any  adaptive  refinement  scheme  is  composed  of  am  error  indicator  amd  a  refinement 
strategy.  Adaptive  remeshing  was  chosen  for  the  refinement  strategy.  This  seems 
natural,  as  it  adlows  to  couple  aidaptation  amd  body  motion  in  a  straughtforwaud  way. 

3.1  Development  of  Suitable  Error  Indicators  for  Adaptive  Remeshing:  We  extended 
the  highly  successful  modified  interpolation  theory  error  indicator  [2,11]  to  3-D.  This 
error  indicator  can  be  generalized  to  multidimensional  situations  by  defining  the  fol¬ 
lowing  tensors: 

/  \Nl^\\Ni\]Uj\dn  ,  (9) 

Ja 

(D*)*!  =  f  \N!l,miUj\dn  ,  =  h?\  f  Nf^NidSl  Uj  \  ,  (10) 

which  yield  an  error  matrix  E  of  the  form: 


E,,  E„ 

E,,  E,,  E,, 


=  X. 


--1 


(11) 


i'x*  Ejfx  Exx 

The  principal  eigenvalues  of  this  matrix  aire  then  used  to  obtaun  reduction  paurauneters 
in  the  three  associated  eigenvector  directions.  Note  that  due  to  the  s3Tnmetry 
of  E  this  is  an  orthogonal  system  of  eigenvectors  that  defines  a  local  coordinate  sys¬ 
tem.  The  required  eigenvalue/eigenvector  solvers  for  3X3  matrices  were  vectorized  and 
incorporated  into  the  code. 

3.2  Smoothing  of  Element  Size.  Stretching  amd  Stretching  Directions:  In  the  context 
of  transient  problems  it  is  very  importamt  to  generate  grids  which  do  not  exhibit  mini¬ 
mum  element  sizes  that  aue  much  smaller  than  the  prescribed  minimum  element  size. 


Practical  calciilatioos  indicated  that  the  grids  produced  by  simply  taking  the  described 
error  indicator  and  the  resulting  distribution  of  element  sizes,  stretchings  and  stretch¬ 
ing  directions  did  not  meet  this  requirement.  In  other  words,  they  were  too  irregular. 
Far  superior  grids  were  obtained  by  smoothing  and  limiting  the  initial  distributions 
obtained  for  element  sizes,  stretchings  and  stretching  directions. 

4.  RIGID  BODY  MOTION 

Rigid  body  motion  algorithms  are  essential  for  moving  body  simulations.  On  the 
other  hand,  just  as  CFD  is  more  than  a  flow  solver,  so  rigid  body  motion  is  more 
than  a  rigid  body  timestepping  scheme.  We  found  that  in  order  to  make  rigid  body 
motion  user-firiendly,  i.e.  accessible  to  a  larger  user  group,  significant  graphical  interface 
developments  were  required. 

4.1  Development  of  a  Rigourous  Rigid  Body  Motion  Algorithm:  If  one  simply  inte¬ 
grates  the  motion  of  rigid  bodies  without  enforcing  the  rigidness  explicitly,  the  shape 
may  change  over  the  course  of  the  calculation.  By  transforming  at  every  timestep  to 
the  body  reference  system,  this  is  avoided.  At  the  seime  time,  the  moments  of  inertia 
tensor  only  needs  to  be  inverted  once,  and  is  much  easier  to  construct. 

4.2  PREMOV;  In  order  to  avoid  costly  mistakes  by  wrong  user-input,  we  developed 
a  graphical  pre-processor  for  cases  where  the  body  motion  is  prescribed.  This  code 
reads  the  same  input  data  as  the  actual  fiow-solver,  but  displays  immediately  the 
trajectories  prescribed.  In  doing  so,  the  user  can  assess  before  the  nm  whether  the 
correct  trajectories  were  specified  in  the  input  file.  Given  that  a  user  may  have  very 
specific  motion-types  in  mind,  a  graphiced  checking  tool  like  PREMOV  is  essential. 

5.  VISUALIZATION  AND  ANIMATION 

Many  of  the  breakthroughs  that  took  place  over  the  performance  period  of  the  present 
grant  were  only  possible  due  to  the  advent  of  powerful  3-D  graphics  workstations.  On 
the  other  hand,  we  had  to  harness  this  power  by  developing  extensive  libraries  to  display 
whatever  was  required.  The  following  tools  were  developed  as  part  of  this  effort. 

5.1  ONLIDISPL:  This  library  allows  to  display  any  desired  quantity  while  the  CFD 
code  is  running  interactively  on  a  workstation.  Thus,  it  elUows  to  monitor  and  check 
a  run  as  it  proceeds  in  real  time.  ONLIDISPL  has  proven  invaluable  for  input  data 
checking  and  debugging,  adding  a  new  dimension  of  user-friendliness.  On  the  other 
hand,  it  is  clear  that  big  runs  can  not  be  done  on  a  workstation.  But  the  user  may 
want  to  nm  a  coarser  mesh  interactively  on  the  workstation,  make  sure  everything  is 
working  as  desired,  and  then  run  the  big  problem  on  a  supercomputer  in  batch  mode. 

5.2  FEPOST3D/MOVIESUBS:  FEPOST3D  is  based  on  FEPLOT4D,  our  standard 
3-D  plotting  capability,  but  performs  all  the  CPU-intensive  filtering  operations  on  the 
supercomputer  within  any  of  the  3-D  codes.  Only  the  plane  or  surfsM:e  triangulations 
jure  sent  back  to  the  SGI-IRIS-4D/80GT  for  plotting.  The  user  specifies  before  the 


run  the  planes  and  surfaces  to  be  inspected.  The  typical  size  of  plot-files  is  reduced 
from  ISO-lSOMbjrtes  (complete  3-D  flowfield,  2Mtetra)  to  2.5-4Mbytes  (plane/surface). 
Although  seemingly  trivial,  this  capability  is  essential  when  trying  to  produce  a  movie, 
or  running  large  3-D  problems  via  a  network. 

5.3  FEMOVIE3D;  We  completed  a  first  movie-making  capablity  for  3-D  runs.  After 
obtaining  the  dumps  from  the  CRAY,  FEMOVIE3D  will  rtm  in  batch  mode  on  the  IRIS 
and  generate  a  series  of  pixel-dumps.  These  pixel-dumps  are  then  shipped  to  a  movie¬ 
making  center  to  obtain  either  VHS  or  0.25in  format  movies.  Although  we  consider 
this  as  only  the  first  generation  of  a  more  sophisticated  movie-making  software,  it 
still  represents  a  significant  development,  as  it  enabled  us  to  assess  the  problems  and 
important  ingredients  that  a  good  movie-making  toolkit  must  possess. 

6.  DEMONSTRATION  CALCULATIONS  AND  RESULTS 

6.1  Unstructured  Grid/Remeshiny  TVansient  .?-T)  Runs:  With  the  developed  tools  we 
performed  the  first  ever  transient  unstructured  3-D  runs  with  adaptive  remeshing.  The 
res^llts  are  summarized  in  [2],  which  is  included  here  as  Appendix  5. 

7.  PUBLICATIONS 

All  of  the  developments  listed  above  were  reported  extensively  in  the  literature.  The 
main  papers  published  are  listed  in  chronologic^d  order: 

[1]  R.  Lohner  and  J.D.  Baum  -  Unstructured  Grid  Methods  for  Store  Separation; 
Proc.  8th  JOCG  Airframe/ Stores  Compatibility  Symp.  ,  Ft.  Walton  Beach,  FI., 
October  23-25,  1990. 

[2]  R.  Lohner  -  Three-Dimensional  Fluid- Structure  Interaction  Using  a  Finite  Element 
Solver  and  Adaptive  Remeshing;  Computer  Systems  in  Engineering  1,  2-4,  257-272 
(1990). 

[3]  R.  Lohner  and  J.D.  Baum  -  Three-Dimensional  Store  Separation  Using  a  Finite 
Element  Solver  and  Adaptive  Remeshing;  .AlAA-91-0602  (1991). 

[4]  J.  Cabello,  R.  Lohner  and  O.P.  Jacquotte  -A  Variational  Method  For  the  Opti¬ 
mization  of  Directionally  Stretched  Elements  Generated  by  the  Advancing  Front 
Method  (AFM)  -  Proc.  Third  Lit.  Conf.  on  Numerical  Grid  Generation  in  CFD 
and  Related  Fields,  Barcelona,  Spain,  June  3-6  (1991). 

[5]  J.  Cabello,  R.  Lohner  and  O.P.  Jacquotte  -  A  Variational  Method  for  The  Opti¬ 
mization  of  Two-  and  Three-Dimensional  Unstructxued  Meshes  -  Proc.  1st  U.S.  Na¬ 
tional  Congress  on  Computational  Mechanics,  Chicago,  Illinois,  July  21-24  (1991). 

[6]  J.  Cabello,  R.  Lohner  and  0-P.  Jacquotte  -  A  Variational  Method  for  the  Op¬ 
timization  of  Two-  amd  Three-Dimensional  Unstructured  Meshes;  AL4A-92-0450 
(1992). 


8.  CONCLUSIONS  AND  OUTLOOK 


We  have  developed  a  CFD  capability  to  compute  3-D  compressible  flows  with  moving 
bodies.  This  is  a  new  capability  that  so  far  is  unmatched.  It  represents  a  new  major 
step  in  CFD,  as  it  allows  to  compute  flows  of  engineering  interest  that  were  previously 
imtractable  by  CFD  means.  The  developments  that  took  place  over  the  last  two  years 
have  been  decribed. 

In  the  future,  we  plan  to  expand  the  developed  capabilities  as  follows: 

a)  Development  of  Flow  Solvers  for  high  Re-number  Viscous  Flows:  We  will  extend 
the  linelet-concept  to  the  fully  coupled  linear  equation  systems  that  arise  for  im¬ 
plicit  Navier-Stokes  solvers. 

b)  Turbulence  Models:  We  will  incorporate  the  k  —  e  model  into  the  implicit  Navier- 
Stokes  solver. 

c)  Development  of  suitable  gridding  algorithms  for  high  Re-number  viscoiis  flows: 
We  plan  to  start  with  a  semi-structured  approach,  whereby  we  construct  prisms 
away  from  the  triangulated  surfaces.  Although  not  general  enough,  this  will  serve 
as  a  starting  point  for  later  developments,  and  a  way  to  gather  experience  for 
simulations  of  high  Re-number  flows. 

d)  Development  of  suitable  error  indicators  for  high  Re-number  viscous  flows:  The 
idea  here  is  to  develop  error  indicators  that  sense  were  to  refine  boundary  layers 
and  shear  layers,  and  that  work  even  for  highly  stretched  grids. 

e)  Improvements  in  movie-making  capabilities:  The  mean  aim  of  this  improvement  is 
to  be  able  to  run  the  movie  on  the  workstation  before  going  to  the  movie-making 
center.  At  the  same  time,  we  must  be  able  to  compress  the  information  to  an 
extent  that  a  3min  movie  can  be  stored  effortlessly  on  a  small  disk.  This  will  be 
accomplished  through  image  compression  algorithms. 

f)  FHirther  test  nins;  we  plan  to  look  for  available  experimental  and/or  numerical 
data  in  the  literature  in  order  to  set  up  runs  to  test  the  algorithms  developed. 
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APPENDIX  1-4:  GRID  OPTIMIZATION 


1  -  Introduction 

Unstructured  meshes,  especially  triangles  in  two  dimensions  and  tetrahedra  in 
three  dimensions,  are  recognized  as  the  more  suitable  to  fit  any  complex  arbitrary 
domain  shape.  Various  algorithms  [1-7,15]  have  been  developed  to  fill  any  arbitrary 
domain  while  preserving  element  quality  to  some  extent.  Unfortunately,  for  most 
techniques,  the  required  versatility  of  the  algorithms  used  to  generate  unstructured 
grids  leads  to  some  loss  of  mesh  quality.  Since  finite  element  solvers  are  sensitive  to 
the  quality  of  the  grid  employed  and  all  generation  techniques  to  date  may  generate 
bad  elements,  means  to  improve  or  smooth  an  initial  mesh  must  be  found. 

Most  of  the  smoothing  techniques  available  in  the  literature  are  based  on  the 
spring  system  analogy  [8,9].  The  mesh  is  viewed  as  a  network  of  nodes  connected 
by  springs  of  constant  stifhiess.  The  smoothing  process  consists  in  moving  the 
nodes  tmtil  the  spring  system  is  in  equilibrium.  Even  though  this  method  is  easy 
to  implement  and  presents  the  advantage  of  requiring  low  CPU-time  and  memory, 
some  drawbacks  have  been  reported.  For  example,  for  non-convex  domains,  the 
spring  analogy  smoothing  may  lead  to  negative  jacobians  in  the  vicinity  of  a  concave 
boundary. 

2  -  Qntimiaation  of  two-dimensional  unstretched  meshes 

In  the  present  work,  we  started  using  a  variational  method  [10-14]  to  optimize 
two-dimensional  xmstretched  triangular  meshes.  The  approach  presented  is  based 
on  the  principles  of  continuum  mechanics  (appendix  1).  The  theoretical  fraihework 
is  inspired  by  the  three-dimensional  non  linear  elasticity  theory.  Starting  from 
basics  geometrical  axioms,  lying  on  mechanical  properties,  we  introduce  a  functional 
defined  as  a  measure  of  the  deformation  between  a  reference  cell  and  a  current 
cell.  This  functional  depends  upon  the  principal  invariants  of  the  deformation 
matrix  (metric  tensor)  derived  from  the  transformation  between  the  two  cells.  The 
mesh  optimization  leads  to  a  minimization  problem  of  em  energy-like  hyperelastic 
material,  quasi  incompressible.  A  convexity  property  is  imposed  to  the  functional; 
then  the  minimization  problem  is  well-posed.  The  minimization  will  be  performed 
by  a  conjugate  gradient  method  (appendix  4). 

3  -  Qntimigation  of  twrydimen-ojonal  Stretched  meshes  (art.  Barcelona  ) 

Dimensional  features  appear  quite  often  in  multi-dimensional  fluid  flows  (e.g. 
shocks,  contact  discontinuities,  etc  ...).  It  hzis  been  shown  by  previous  work  [7]  that 
allowing  directionally  stretched  elements  is  an  efficient  way  to  increase  the  solution 
accuracy  without  drastically  increasing  the  number  of  elements. 

The  advancing  front  method  (AFM)  proposed  by  Permre  et  al.  [7],  Lohner  [15] 
,  Lohner  and  Parikh  [6]  is  able  to  generate  stretched  elements  in  the  direction  of  the 
flow  features.  Unfortimately,  the  unstructured  grid  generated  with  the  advancing 
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front  technique  may  lack  the  desired  distribution  of  the  grid  points  in  some  regions. 
As  the  ntunerical  stability  and  accuracy  of  flow  solutions  may  be  affected  if  the  grid 
points  are  not  distributed  smoothly,  the  imtial  mesh  generated  by  the  AFM  must 
be  smoothed. 

The  variational  method  was  extended  to  the  optimization  of  two-dimensional 
stretched  triangular  meshes  generated  by  the  AFM.  Being  able  to  cope  with  the 
optimization  of  two-dimensional  unstretched  meshes,  we  needed  to  devise  a  way 
to  transform  the  stretched  mesh  into  an  unstretched  mesh.  This  was  done  by  the 
introduction  of  a  transformed  space.  A  transformation,  depending  on  the  mesh 
parameters  used  to  construct  a  stretched  element,  was  performed  to  obtain  an 
element  which  shotild  look  like  an  equilateral  element.  This  transformation  is  just 
a  rotation  of  the  element,  in  order  to  align  the  stretching  direction  with  the  x-axis 
and  a  shrinking  of  the  element  (along  the  x-axis)  of  a  factor  equal  to  the  inverse 
of  the  stretching  factor.  Asuming  that  the  element  generated  in  the  physical  space 
has  exactly  the  characteristics  given  by  the  mesh  parameters,  then  the  resulting 
element  after  transformation  would  be  an  equilateral  triangle.  The  optimization  is 
performed  on  the  transformed  space.  Finally,  the  optimized  mesh  is  transformed 
back  (using  the  inverse  transformation)  from  the  transformed  into  the  physical 
space. 

4  -  Optimiigation  of  three-dimensional  stretched  meshes 

In  three  dimensions,  a  transformation  between  the  physical  and  transformed 
space  was  devised  as  an  extension  of  the  two-dimensional  transformation.  'In  the 
transformed  space,  all  tetrahedra  are  assumed  to  be  regulars  (all  faces  being  equi¬ 
lateral  triangles).  The  variational  method  is  applied  in  the  transformed  space. 
First,  a  measure  of  the  deformation  between  a  transformed  and  reference  element  is 
computed.  The  total  deformation  is  obtained  by  stimmation  of  the  elementary  con¬ 
tributions.  Second,  the  interior  mesh  points  are  repositioned  in  order  to  minimize 
the  total  deformation.  Finally,  the  elements  are  transformed  back  to  the  physical 
space. 

5  -  Optimi7.«.t.ion  of  three-dimensional  iinstretched  meshes  (art.  AIAA-92-0450  ) 

For  the  special  case  of  three-dimensional  unstretched  tetrahedral  meshes  there 
is  no  need  to  go  back  and  forth  between  the  physical  space  and  the  transformed 
space.  The  method  is  directly  applied  in  the  physical  space.  A  special  version  of 
the  code  has  been  written  for  the  optimization  of  unstretched  meshes.  This  version 
requires  less  memOTy  and  less  cpu-time  than  the  general  one. 

6  -  Post-processing 

Using  a  moving  node  method  we  are  tied  to  the  topology  given  by  the  initial 
mesh.  Some  routines  have  been  written  which  change  the  topology  of  the  mesh. 
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These  routines  are  used  to  either  remove  cells  with  negative  jacobians  or  cells  which 
are  marked  as  distorted  according  to  a  quality  criterion.  These  routines  have  been 
isolated  and  can  be  used  as  a  “black  box”  during  the  post-processing  step  of  any 
initial  or  optimized  mesh.  During  the  conjugate  gradient  iterations  we  may  get 
folded  elements  with  negative  jacobians.  These  elements  will  be  unfolded  as  the 
conjugate  gradient  iterations  are  carried  along  and  no  cell  with  negative  jacobian 
will  remain  when  the  minimum  is  reached.  Nevertheless,  these  negative  jacobians 
will  usually  appear  in  some  pathological  regions  of  the  mesh.  We  suspect  that 
these  elements  originate  where  the  advancing  front  closes  small  gaps  in  order  to  fill 
the  whole  domain.  Therefore,  it  is  preferable  to  take  out  these  elements,  chamg- 
ing  locally  the  topology  of  the  mesh,  rather  than  carrying  the  conjugate  gradient 
iterations  to  convergence. 

7  -  MiniTni7.ation  algorithm 

Within  the  domain,  the  optimal  location  of  the  physical  nodes  is  computed 
by  an  iterative  conjugate  gradient  algorithm  (appendix  4).  The  multidimensional 
minimization  problem  is  reduced  to  a  succession  of  one  dimensional  problems  in  a 
descent  direction.  The  one  dimensional  search  for  a  minimmn  leads  to  find  the  root 
of  a  three  (in  two  dimensions)  or  five  (in  three  dimensions)  degree  polynomial  of 
one  ''?jiable.  The  coefficients  of  this  polynomial  depend  on  the  differential  of  cr. 
The  formtilas  for  these  differentials  are  given  in  appendices  2  and  3. 

8  -  Conclusion 

The  scope  of  a  variational  method,  widely  used  for  the  optimization  of  struc¬ 
tured  grids,  has  been  extended  in  order  to  perform  the  optimization  of  two-  and 
three-dimensional  stretched  or  unstretched  imstructured  meshes.  This  moving  node 
method  is  not  tied  to  the  advancing  front  method  or  to  any  pzurticular  generation 
technique  and  can  be  initialized  by  any  arbitreury  grid,  even  a  grid  with  overlapped 
cells.. 

In  the  paper  presented  in  Barcelona,  the  variational  method  was  able  to  give 
acceptable  results  in  a  two-dimensional  case  where  the  spring  analogy  failed.  The 
robustness  and  efficiency  of  the  method  was  tested  for  a  mesh  in  which  the  location 
of  the  points  was  altered.  The  method  was  able  to  unravel  this  “chaotic”  mesh 
and  to  give  a  smooth  one.  For  stretched  meshes  the  method  was  able  to  exert 
more  control  aa  the  grid  according  to  the  parameters  specified  by  the  user  in  the 
background  mesh. 

In  the  paper  presented  at  the  AIAA  meeting,  for  the  optimization  of  un¬ 
stretched  three-dimensional  tetrahedral  meshes,  a  more  appropriate  formulation 
than  the  general  one,  have  been  presented  in  order  to  save  cpu-time  and  memory 
space.  For  all  the  examples  presented  the  quality  of  initial  meshes,  already  smoothed 
by  the  spring  analogy,  have  been  improved  with  an  increase  of  about  eight  to  nine 
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per  cent  around  the  optimal  angle  (angles  between  fifty  and  eighty  degrees)  along 
with  a  decrease  of  angles  leading  to  fiat  elements  (i.e.  angles  lower  than  thirty  and 
greater  than  one  htmdred  and  ^y  degrees).  Also,  for  all  examples,  an  increase  of 
a  factor  two  approximately  is  obtained  for  the  minimum  radii  ratio  of  the  inscribed 
over  the  circumscribed  sphere.  This  is  an  important  result,  given  that  the  minimum 
radii  ratio  drives  the  allowable  time  step  of  finite  element  solvers. 
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APPENDIX  1 


Construction  of  the  functional 
Theory 

We  review  the  mechanical  and  mathematical  properties  of  the  functional  used  to 
measxire  the  deformation.  A  more  detailed  explanation  of  the  foundations  of  the 
method  can  be  found  elsewhere  [10,14]. 

1.1.  -  Mechanical  Assumptions 

Let  us  consider  a  transformation  x(£)  between  a  reference  element  and  a  physical 
element.  First,  we  restrict  the  dependence  of  the  functional  up  to  the  first  order 
derivative  of  the  transformation  : 

<r  =  <T(x,F)  with  F  =  Vx  (1) 

As  the  functional  should  not  change  if  we  apply  a  translation  to  the  physical  element, 
we  further  restrict  a  to  rely  on  the  transformation  x  solely  through  its  gradient  F. 

a(x,F)  =  <T(F)  (2) 

Since  any  rotation  of  the  physical  element  must  not  change  the  measure  we  require 
the  following  identity  to  hold  for  <t  :  , 

o(F)  =  <t(Q  .  F)  for  any  orthogonsJ  matrix  Q  (3) 

Furthermore,  since  the  measure  must  also  be  invariamt  for  any  reference  rotation 
we  impose  : 


<r(F)  =  <7(F  •  Q)  for  any  orthogoneJ  matrix  Q  (4) 

To  sum  up,  we  reqviire  the  functional’s  invariance  to  rigid  motions  of  both  reference 
and  physical  element.  It  means  that  a  functional  following  the  four  axioms  and 
properties  aforementioned  is  independent  of  the  orthonoimal  basis  in  which  the 
gradient  F  is  computed.  This  property  can  be  written  as, 

(T  =  <T(q.F.Q‘)  for  any  orthogoneJ  matrixQ  (5) 

The  axioms  and  properties  (1)  to  (4)  zue  well  known  in  three-dimensional  elasticity 
theory  [lOj.  They  characterize  the  energy  of  an  hyperelastic  (1),  isotropic  (4). 
homogeneous  (2)  material  satisfying  the  zuciom  of  frame  indifference  (3). 

From  these  properties,  it  can  be  shown  that  the  functional  a  depends  only  on 
the  invariants  of  a  matrix  C,  called  the  right  Cauchy-Green  tensor  of  the  trsinsfor- 
mation  x,  and  defined  as 


c  =  •  F  =  Vx*  •  Vx 


(6) 


At  zmy  point  C  =  (^>»7>C)  in  tiie  reference  space,  the  invariants,  noted  Ii,l2  and  I3 
are  given  by  the  formulae  : 


11  =  tr  (C) 

=  l|F||^ 

=  l|x<lll  +  l|xvllv  +  llx<||^  (7) 

12  =  tr(Cof(C)) 

=  ||CofF||^ 

=  ||x^  A  x,||^  +  ||x<  A  X<:ll2  +  ||X,  A  x^Wl  (8) 

13  =  det  (£) 

=  (detF)^ 

=  J  ’  ^  (9) 

with. 


J  =detF 

=  det(X(,x„Xc)  (Jacobian)  (10) 

and  the  symbols  ll.||m»  IMIv*  tr  (.),  det  (.)  and  Cof(.)  denoting,  respectively,  the 
matrix  norm,  the  vector  norm,  the  linear  operator  trace  and  the  two  non  linear 
operators  determinant  and  cofactor. 

<7  =  <7(Il,Ij,l3)  (11) 

In  order  to  oootiol  the  element  orientation  and  prevent  negative  j2u:obians,  the 
invariant  I)  (insensitive  to  the  orientation)  is  replaced  by  the  Jacobian.  From  now 
on,  we  are  seeking  for  a  functional  in  the  general  form. 


<T  =  <7(Ii,l2,  J  ) 


(12) 
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Vowing  to  establish  a  well-posed  minimization  problem,  we  force  the  functional  to 
be  locally  convex  in  the  neighborhood  of  a  rigid  transformation. 

1.2.  -  Mathematical  Assumptions 

Rigid  treinsformations  -  i.e.  translations,  rotations  -  verify  the  following  equivalent 
properties  : 


x(^)is  a  rigid  transformation  (13) 

F  =  Vx  is  a  direct  orthogonal  matrix  (14) 

F=  CofF  and  J  = -f-1  (15) 

(Ii,l2,  J)  =  (3,3,l)  (16) 


In  the  sequel,  rigid  transformations  will  be  characterized  by  the  subscript  index  0. 

Owing  to  the  shape  preserving  property  of  a  rigid  transformation  we  impose  a 
normalization  condition  : 


<t|o  =  0  (17) 

Then,  an  equilibritim  condition  is  settled  for  rigid  transformations, 

D<t|o  =  0  (the  null  tensor)  (18) 

Finally,  in  an  attempt  to  obtain  a  unique  minimum,  the  functional  is  <issumed  to 
be  convex  in  the  neighborhood  of  a  rigid  transformation  ; 

D^<t(o  >  0  (positive-definite  tensor)  (19) 

In  these  expressions  D  and  denote  respectively  the  first  and  second-order  differ¬ 
ential  of  the  functional. 

1.3.  -  Functional  Chosen 

We  are  now  in  position  to  move  on  to  the  practical  choice  of  the  functional.  Several 
candidates  verifying  both  mechanical  and  mathematical  requirements  can  be  found 
in  the  literature  [11].  In  order  to  devise  a  functional  more  amenable  to  numerical 
computations  we  choose  a  polynomial  form,  discussed  in  details  in  [5],  that  can  be 
expressed  in  three  dimensions  as  : 


<T34  =!  C  i(Ii  —  I3  —  2)  -|-  C  2(12  —  2I3  —  1)  -f-  C  3(  J  —  1)^ 


(20) 
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with  the  convexity  condition  (19)  leading  to  the  inequality 

3C3>4(Ci  +  C2)>0  (21) 

In  two  dimensions,  the  expression  reduces  to 


9 


APPENDIX  2 


Two-Dimensional  Optimization 
Numerical  aspects 

The  numerical  approximation  of  the  trsmsformation  and  its  derivatives  are 
presented. 

introduction 

We  w£mt  to  compute  the  deformation  between  a  current  element  in  the  physiczd 
space  (  X,  y  )  and  a  reference  element  in  a  reference  space  (f,  r)).  The  current 
element  is  a  triangle  given  by  its  vertices  coordinates  ((xi);  i  =  1,3)  while  the 
reference  element  is  a  triangle  given  by  its  vertices  coordinates  ((^.);  i  =  1,3).  A 
linear  transformation  x  is  assumed  between  the  reference  triangle  and  the  cmrent 
triangle. 

2.1.  -  Linear  transformation  x  : 

For  the  vector  transformation  x  we  assiune  that  each  sczdar  transformation,  i.e.  i 
and  y,  is  a  linear  function  of  ^  r;)  i.e.  : 


x(i)  =  W  6  »?),  y(  ^  T?))‘ 

X  =  +  brj  +  c  a  =  (a*,  »  b  =  (6*,  6^)  ,  c  =  (cx,  Cy) 

Following  the  line  of  the  classical  finite  element  method  the  above  vector  transfor¬ 
mation  can  also  be  defined  by  : 


i=3 

X  =  ^  (2) 

i=t 

Knowing  the  vertices  coordinates  in  the  physical  space  the  vector  transformation  is 
completely  defined  as  soon  as  the  element  shape  functions  iV,  ,  i  =  1,3  are  known. 

2.2.  -  Element  Shape  Functions  : 

From  expression  (2),  the  linear  shape  functions  are  defined  such  that  Ni  has  the 
value  iinity  at  node  i  and  is  zero  for  the  remaining  nodes  : 


^i(  (j  V)  =  +  fti'?  +  Ci  1  <  *  <  3 


(3) 

(4) 
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with  6ij  referring  to  the  kronecker  symbol  and  Oi,  bi  and  Ci  being  scalar  values. 
2.3.  -  Computation  of  the  Shane  Functions  : 


From  equations  (3)  and  (4)  we  deduce  that  the  tinknowns  (a,,  hi,  Ci)  axe  solutions 
of  the  system  : 


'6  »?i  1' 

Oi 

'6ii' 

6  m  1 

hi 

= 

Si2 

.6  »73  1 , 

^Ci 

hi3 

That  we  cam  rewrite, 


(5) 


with, 


A  •  u  =  r 


(6) 


'6  m  i' 

cii 

'Sii' 

A  = 

6  1 

;  u  = 

b. 

;  r  = 

6i2 

.6  1 . 

_Ci 

The  system  of  equations  defined  by  (6)  has  a  unique  solution  if  the  determinant  of 
the  matrix  A  does  not  vanish.  Introducing  the  circular  permutation  (ij,k)  between 
indices  1,2,3  we  obtain  : 


Oi  =(»?/ -»?*)/ deter 

=((k-^j)/ deter  (7) 

Ci  =((jrjk-Unj)/  deter 


with 


deter  =  det  A 

=  2(  area  of  the  reference  element  ) 


(8) 


Remark  2  : 

-  The  area  of  the  reference  triangle  is  defined  by 

area  =  ^  |  (6  -  ^  (6  “  6)1 

where  the  symbol  A  represents  the  vector  product. 

From  expression  (4)  amd  (7)  the  element  shape  functions  aure  completely  defined. 
Afterwau'ds,  the  vector  transformation  is  defined  by  expression  (2). 
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-  One  convenient  and  practical  way  to  verify  that  the  basis  functions  are  correctly 
defined  consists  in  noticing  that  the  stim  of  the  three  basis  functions  is  a  linear 
function  taking  the  value  unity  on  three  points.  Consequently,  the  siun  Ni  + 
N2  +  N3  is  the  constant  function  unity  and  we  have  : 

Ni-\-N2  +  N3  =  1  (9) 


We  will  need  to  compute  the  derivatives  of  the  vector  function  x.  According  to 
expression  (2)  the  derivatives  are  : 


dNi  dN2  dNi 

Xt  =  X.-g^+X,-g^+X3-g^ 
»  4-V 


(10) 


2.4.  -  Computation  of  the  derivatives 

The  elements  shape  functions  being  linear,  their  derivatives  are  constants.  These 
constants  are  given  by  derivation  of  formula  (4)  leading  to  : 


dNi 


=  Ot 


dNi 


=  6i 


(11) 


Remark  4  : 


-  Using  expression  (9)  we  deduce  that  the  sum  of  all  the  derivatives  vanishes 
'  dNi  ,  dN2  ,  dN3  dNi  ,  dN2  ,  dN3  ^ 

■^  +  'ar  +  ‘ar  "  "aS"  aT  “  *  ’ 

The  equalities  (12)  offers  a  good  way  to  verify  numerically  the  computation  of 
the  element  shape  function  derivatives  . 

At  this  point,  we  defined  all  the  elements  required  to  compute  the  transformation 
amd  its  derivativea.  The  forthcoming  section  is  devoted  to  the  computation  of  the 
invariants  for  the  right  Cauchy-Green  tensor  and  their  derivatives  which  will  be 
used  for  the  minimization  of  the  functionsil.  The  computation  of  the  functional  and 
its  derivatives  will  be  used  during  the  minimization  performed  by  the  conjugate 
gradient. 

2.5.  -  Computation  of  the  invariants 

Let  us  remind  that  we  have  chosen  to  approximate  the  vector  transformation  by  a 
linear  function  : 
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with 


X  =  xi  JVi  +  X2N2  +  X3N3 


the  derivatives  are  constants  given  by  the  expressions  : 


r  x^(  X  )  =  xifli  +  X2a2  +  xafls 
\  x,(  X  )  =  X161  +X262  +X363 


(13) 


Remark  5  : 

-  FVom  expression  (13),  the  derivatives  are  linear  function  of  the  vertices  coordi¬ 
nates  Xi,X2,X3 


The  vector  function  gradient,  F,  is  a  (  2  x  2  )  matrix  defined  by 


F  =  Vx 


=  [xc  .  X,) 

(14) 

The  Cauchy-Green  matrix  is  : 

C  =  F‘  F 

(15) 

Finally,  the  two  invariants  of  the  Cauchy-Green  tensor  are  : 

h  =  tr  [C] 

(16) 

I3  =  det  (C)  =  J  2 

with 

J  =  det  (F)  , 

(17) 

the  Jacobian  of  the  transformation 
2.5  -  Computation  of  the  invariants 


2.5.1.  -  Computation  of  the  iacobian 
The  Jacobian  is  computed  xising  the  expression  : 

J  =  det  (X(,  x,)  (18) 


Remark  6  : 
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1  -  The  Jacobian  is  independent  of  ^  and  17.  It  only  depends  on  the  vertices  coor¬ 

dinates  of  the  current  element  (physical  space). 

2  -  The  Jacobian  is  an  homogeneous  polynomial  of  degree  2  i.e. 

J  =  aijXiyj  with  Oij  e  ^  (19) 

l<i<3 

1<J<3 


where  Xi,yi  represents  the  cartesian  components  of  vertex  Xj. 

3  -  The  Jacobian  represents  the  ratio  of  the  current  element  surface  over  the  refer¬ 
ence  element  surface. 

2.5.2.  -  Computation  of  the  first  invariant 
By  definition,  the  first  invariant  is  : 

h  =  tr  [C]  (23) 

thus,  we  deduce  from  expressions  (14)  and  (15), 

h  =  -h  ||x,||^  (24) 

which  allows  us  to  compute  Ii  through  expression  (13)  of  the  vector  function  deriva¬ 
tives. 

Remark  7  ; 

1  -  The  first  invariant  is  independent  of  ^  auid  77,  depending  only  on  the  vertices 

coordinates. 

2  -  Ii  is  a  second  order  homogeneous  polynomied  of  the  coordinates  (xi ,  j/i ,  12 ,  y? , 

X3,  ya)- 

2.6.  -  Differentiation  of  the  invariants 

So  far,  we  have  seen  how  the  two  invariants  are  computed.  From  expressions 
(18)  and  (21),  the  computation  of  the  invariants  J  and  Ii  is  related  to  the  com¬ 
putation  of  the  vector  functions  X(  and  x,,.  Now,  we  are  going  to  give  the  formulas 
for  the  differentiation  of  the  invariants  and  their  derivatives. 

2.6.1.  -  Notations  : 

In  the  following  we  note 

X  =  (xi,  yi,  X2,  y2,  ®3,  ys)  =  (Xi*,  X2*,  X3*)  (22) 


and 
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b  =  (hi  hi  hi  hi  h’)  =  (h  h  2*,  b  ,')  (23) 

with, 


X,  = 


Xi 

Vi 


and 


representing  respectively  the  vertices  of  the  triangle  and  a  direction  vector  at  the 
vertices  (This  direction  vector  may  be  the  descent  direction  at  the  vertex). 


Remark  8  : 

-  The  3li*-vector  x  is  the  variable  of  our  problem  (x  represents  the  three  vertices 
coordinates).  According  to  expression  (2),  the  vector  tremsformation  x  is  com¬ 
pletely  defined  when  the  three  vertices  coordinates  Xi ,  X2  and  X3  are  defined. 
Thus,  there  is  a  one-to-one  correspondence  between  the  vector  transformation 
X  and  the  S*  vector  x. 


2.6.2.  -  Differentiation  of  th«»  transformation  derivatives 
From  expression  (13)  we  get  for  the  vector  function  derivatives  : 


■X{(X1,X2,X3) 

Xij(2;i,X2,X3) 

.  y{(yi,y2,y3) . 

X,  — 

.  y,(yi,y2,y3) . 

The  above  expressions  show  the  derivatives  dependence  on  the  variable  x.  Using 
the  standard  dififerentiation  rules,  the  vector  function  differential  is  : 


and. 


with. 


dx((x)-  h 


d  I^(xi,X2,X3)  •  (ftf,/if,^f) 

.  d  y«(yi,y2.y3)  • 


dx,(x)-  h 


d  x,(xi,X2,X3)  • 

d  y^Cyi.ya.ys)  •  {h\,h^,hl) 


(24) 


(25) 


d  Xj(*l,X2,X3)  •  (/lf,/»2,/lf)  =  ^^^3 


(26) 


and  similar  formulas  holding  when  interchanging  x  by  y  and  the  subscript  ^  by  77. 


2.6.2.  -  Differentiation  of  the  iacobian 

From  definition  (18)  of  the  Jacobian  and  using  the  differentiation  rules  we  obtain 
d  J  (x)  •  (  h  )  =  det  (  d  X{(x)(  h  ),x,)  -i-  det  (X{,  d  x,(x)(  h  ))  (27) 
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Remark  9 

1  -  As  the  vector  functions  X(  and  are  linear  functions  of  the  variable  x  we  have 
the  following  relation  : 


d  x^(x)  •  (  h  )  =  x«(  h  ) 

=  h^  (28) 

where  h  ^  is  the  vector  function  defined  by  equation  (13)  when  replacing  x  by 
h  .  Then,  we  can  rewrite  (27)  in  the  form  : 

d  J  (x)  •  (  h  )  =  det  (  h  (,x,)  +  det  (x^,  h  ,)  (29) 

2  -  Using  the  differentiation  rule  we  get  the  higher  order  differential  of  the  jacobiein. 
2.6.4.  -  Differentiation  of  Ii 

Using  the  same  notations  and  applying  the  differentiation  niles  to  formula  (21)  we 
deduce  : 


d  Ii(x)(  h  )  =  2(  h  ^  .  X(  +  h  ,  •  X,)  (30) 

2.2.3.  -  Computation  of  the  functional 

<r  ^  Ci(Ii  -  I3  -  1)  +  C2(  J  -  1)2  (31) 

-  Ii  is  a  polynomial  of  degree  2. 

—  J  is  a  polynomial  of  degree  2. 

—  I3  is  a  polynomial  of  degree  4. 

From  all  the  previous  sections  we  are  able  to  compute  numerically  a  zind  its  deriva¬ 
tives. 
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APPENDIX  3 


Three-Dimensional  Optimization 
Numerical  aspects 

The  numerical  approximation  of  the  transformation  emd  its  derivatives  are 
presented  in  three  dimensions. 

introduction 

By  the  same  token,  following  the  same  line  of  thinking,  the  numerical  approximation 
presented  in  two  dimensions  can  be  extended  to  three  dimensions.  We  waint  to 
compute  the  deformation  between  a  current  element  in  the  physical  spsure  (  x,y,  z 
)  and  a  reference  element  in  a  reference  sp2w:e  (^,  77,  C)-  The  current  element  is 
a  tetrahedron  given  by  its  four  vertices  coordinates  ((  ^  ,);  i  =  1,4)  while  the 
reference  element  is  a  tetrahedron  given  by  its  vertices  coordinates  *  =  1»3)- 

A  linear  transformation  x  is  zissumed  between  the  reference  triangle  and  the  current 
triangle. 

3.1.  -  Element  Shane  Functions  : 

The  basis  functions  Ni  are  linear  functions  defined  by  the  relations  : 

-  ^ij  V»,j;  l<t,j<4  '  (1) 


Ni(i,rj,0  =  ai^  +  biT}  +  CiC-i-di  (2) 

with  =  iijiVjiCj)  til®  coordinates  of  the  vertices  of  the  tetrahedron  in  the  refer¬ 
ence  space  (^,»7,C)* 


3.2.  -  computation  of  Ni(E,ii.  CTl 

The  unknowns  aj,fri,Ci  and  di  are  solutions  of  the  linear  system  : 


{<*1^1  +  4-  d,  =  6ii 

+  biri2  +  c.Cj  +  dj  =  6i2 
<»i^3  +  +  c.Cs  di  =  5i3 

+  CtC4  +  d,  =  ^i4 

The  determinant  of  this  system  is  given  by  : 


deter  = 


6 

'71 

Cl 

Sil 

C2 

V2 

C2 

6i2 

'73 

C3 

^i3 

^4 

'74 

C4 
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If  deter  ^  0  then  the  system  S  has  a  unique  solution. 


Remark  1  : 

-  The  summation  of  the  basis  functions  is  a  linear  polynomial  of  degree  1,  equal 
to  1  at  four  points.  Therefore,  this  polynomial  is  the  constant  polynomial  equal 
to  1  and  we  have  : 


Ni+N2+N3  +  N^  =  1 

3.3.  -  Computation  of  the  derivatives 

All  derivatives  are  constant  and  given  by  the  expressions  : 


dNi 

dNi 

=  h 

dNi 

=  ai 

’  dri 

’  dC 

=  Cl 

dNi 

dNi 

dNi 

=  02 

’  dr, 

=  h 

’  dC 

=  Ci 

dNz 

dNi 

=  &3 

dNi 

di 

-  03 

’  dr, 

’  <9C 

=  C3 

dN4 

dN4 

dN4 

=  04 

’  dr. 

=  h 

’  ac 

=  C4 

Remark 

-  The  derivatives  of  the  sum  of  the  basis  functions  is  null  (i.e.) 

Ea£f=»  fori  =  l,2,3  with({■,^^^“)  =  (f,,,C) 

>=1  ' 

3.4.  -  Computation  of  the  invariants 

The  transformatiom  x  may  be  rewritten  as  : 

1=4 

x  =  x,N,  (3) 

i=l 

with, 


X 


;«  =  1,2, 3,4 


The  derivatives  are  constants  given  by  the  expression  : 
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fori  =  1,2,3  with  =  (C,7?,C)  (4) 


tsl 


wd. 


dNi 

de> 


representing  the  derivative  of  the  shape  function  Ni  with  respect  to  the  variable  . 
We  obtain, 


X((  X  )  =  OiXi  +  02X2  +  03X3  +  04X4 

x,(  X  )  =  5iXi  +  62X2  +  ^3X3  +  64X4  (5) 

x<(  X  )  =  ciXi  4-  C2X2  +  C3X3  +  C4X4 

The  vector  function  gradient,  P,  is  a  (3  x  3)  matrix  defined  by 

F  =  Vx 
=  [x^  ,  x„  xj 

The  Cauchy-Green  matrix  is  :  » 

C  =  F‘  F  (7) 

Finally,  the  three  invariants  of  the  Cauchy-Green  tensor  are  : 


Ii  =  tr  (£) 

=  m\}n 

=  l|3'<llv  +  l(x*»llv  +  ll*<llv  (8) 

la  =  tr(Co£(£)) 

=  l|CofF||^ 

=  ||x<  A  x,,i|^ -H  llx^  A  X<||J-t-||x,  A  x^ll^  (9) 

I3  =  det  (£) 

=  (detF)^ 

=  J  ^  (10) 
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with, 


J  =  detF 

=  det(X{,x„X()  (Jacobian)  (11) 

and  the  symbols  ||.|[m)  II-||v>  tr  (.),  det  (.)  and  Cof(.)  denoting,  respectively,  the 
matrix  norm,  the  vector  norm,  the  linear  operator  trace  and  the  two  non  linear 
operators  determinant  and  cofactor. 

3.5.  -  Differentiation  of  the  invariants 

The  computations  of  the  invariants  are  carried  out  using  formulas  (3)  of  the  basis 
functions  derivatives.  The  next  step  is  to  find  the  formulas  for  the  differential  of 
these  invariant  which  will  be  used  to  compute  the  differential  of  the  functional. 

3.5.1  -  Notations ; 

In  the  following  we  note, 

X  =  (xi,X2,X3,a:4,yi,y2,y3,y4,2l,«2,«3,«4) 


and 


h  ==(ht,hlhlhlh\,hlhlhlhl,hlhlhl) 


with, 


X  i  = 

Xj' 

Vi 

h.= 

\h!] 

[Ml 

representing  respectively  the  vertices  of  the  tetrahedron  and  a  direction  vector  at 
these  vertices. 

3.5.2  -  Differentiation  of  the  tranaformation  derivatives  : 

The  derivative  of  the  transformation  x  is  : 


X^(xi , X2 , X3, I4 ) 

y{(yi,y2,y3,y4) 

X((^i, 22,^3,24)  . 


For  the  differential  of  the  vector  function  x^  at  point  x  applied  to  vector  h  we 
have  ; 
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dx5(x)-  b 


d  i{(x)  •  h 
dy^Cx)-  b 
d  Z{(x)  •  h 


aad  the  same  relations  holding  for  r}  and  C- 
3.5.3  -  Differentiation  of  the  iamhian 
The  Jacobian  being  related  to  the  derivatives  by  the  expression  : 


J  (X)  =  (X{,X„X,;) 

with  X{,Xq,X(,  linear  functions  of  the  variable  x,  the  differentiation  gives  ; 


d  J  (  X  )  •  (  h  )  =  d  X((  X  )(  h  )  •  (x,  A  xc)+ 

d  x,(  X  )(  h  )  •  (xc  A  xe)+ 

d  x<(  X  )(  h  )  •  (X{  A  X,)  (12) 

The  scalar  obtained  above  is  the  stunmation  of  the  scalar  product  between  the 

differential  of  the  derivative  function  and  the  vector  product  of  the  two  remaining 
derivatives.  The  differentiation  of  the  derivatives  simplifies  to  : 

d  X{(  X  )  •  (  h  )  =  x^{  h  ) 

=  h^  (13) 

where  b  ^  is  the  vector  defined  by  relation  (3)  while  replacing  variable  x  by  h  . 
The  differential  of  the  Jacobian  rewrites  : 

d  J(x)-(b)  =  b<-(x,  A  x^)+  b,  (xc  A  X()+  bc-(^€  ^  X,)  (14) 

We  can  also  deduce  the  higher  orders  differentiation  of  the  functions^  J  and  with 
the  same  notations  and  considerations  introduced  previously  we  obtain  : 


d’j(  X  )  •  (  b  ,  b  )  =2(  (X(,  b  ,,  b  c)+ 

(  b  e,x,,  b  c)+ 

(b(,  b,,xc)  I  (15) 
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<1  ’  J  (  X  )•(  b  ,  b  ,  b  )  =  6(  b  b  ,,  b  c)  (16) 

3.5.4.  -  Differentiation  of  L 

Ii=xJ+xJ  +  x2  (17) 

<ili(x)(b)  =  2  (  b  (.X(  +  h  ,.x,  +  h  ^.xc)  (18) 

d^Ii(x)(b,  b)  =  2  (  b(.  bf  +  b  ,.  b  ,  +  b  c- b  c)  (19) 

3.5.5.  -  Differentiation  of  lo 

I2=(X(  A  x,)^+(x,  A  X()2+(x^  A  X{)2  (20) 


d  I2  =  2  (  (X(  A  X,)  •  {hj  A  X,  +  Xj  A  h^} 

+  (x,  A  x<)-{b,  A  x^+x,  A  h^} 

+  (x<  ^  x«)  {b<;  A  x^+x^  A  h^}  )  (21) 


d  ^la  =  2  (  (b(  A  X,  +  x^  A  h,)2 

+  (b,  A  X(.+X,  A  h^)^ 

+  (b<  A  x^  +  x<.  A  hf )^ 

+2  [  (x^  A  x,)  (h^  A  h,) 

(x,  A  x^)  (h,  A  h^) 

(x<  A  x^)-(b^  A  h()  ]  )  (22) 


d  *I2  =  12  (  ((b^  A  X,  +  x^  A  h,)  •  (h^  A  h,) 

+  (b,  A  X^+X,  A  h<.)  (h,  A  b^) 

+  (b<  A  x^  +  x<.  A  h()  (h^  A  b{)  ]  )  (23) 

From  expression  (23)  we  see  that  the  scalar  function  d^l2  is  linear  and  then  for 
the  computation  of  d  (constant  i.e.  independent  of  the  variable  x  )  we  use  the 
previous  expression  replacing  x  by  h  . 
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3.5.6.  -  Differentiation  of  the  functional 


o  =  Ci(Ii  —  I3  —  2)  +  ^2(12  2I3  —  1)  +  C3(  J  —  1)^  (24) 

•  Ii  is  a  polynomizd  of  degree  2  (with  respect  to  the  variable  x  ) 

•  I2  is  a  polynomial  of  degree  4  (with  respect  to  the  variable  x  ) 

•  J  is  a  polynomial  of  degree  3  (with  respect  to  the  variable  x  ) 

Therefore  the  functional  <7  is  a  polynomial  of  degree  6  of  the  vertices  coordinates. 
Reordering  and  gathering  the  terms,  a  rewrites  in  the  gener2d  form  : 

<7  =  C*!!!  +  C2I2  d"  T  J  ^  —  2C3  J  -f  ^ 

with, 

'Y  =  C3  —  Cj  —  2C2  and  h  —  C3  —  2Ci  —  C2 

The  differentiation  of  a  gives  : 

d<7(  X  ).  b  =C\  d  Ii(x)-  h  +C2d  l2(x)  -  h  + 

(27  J(x)-2C3)dJ(x)  -  h 

d  V(  X  ).(  b  ,  h  )=  Cl  d  I  i(  b  )•  b 

+  C2  d  /2(  X  )  •  (  h  ,  h  )  +  27(dJ(  X  )  •  b  )^+ 

[27  J  (  X  )  -  2C3]  d  J  (  X  )  •  (  b  ,  b  ) 

d  V(  X )  •  (  b  ’)  =C2  d  1 2(  X  )  •  (  b  ,  b  ,  b  )+ 

47(dJ(x).  h  )•(  d  J(x)-(h,  h))+ 

[27  J  (  X  )  -  2C3]  d  J  (  b  )  •  (  b  ,  b  )+ 

27dJ(  x)  -  b  d  J(x)-(b,b) 

U  V(«Mb")=C2d  l2(b)(b,  b,  h)+ 

87(  dJ(  X  )  •  h  )  •  (  d  J  (  h  )  •  (  b  ,  b  )  )+ 

67(  d  J  (  X  )  •  (  h  ,  h  ))2 

dV(x).(b®)  =  "'i7d  j(x)-(h,  h)d  J(b)-(b,  b) 
d  *<7(  X  ).(  b  ®)  =  207  d  J  (  h  )  •  (  h  ,  b 
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appendix  4 


Coigugate  Gradient  Algorithm 

V  = 

nunk  =  ndiinn^npoin  (unstretched  meshes) 

=  ndimn*nnodes*nelem  (stretched  meshes) 

<r  is  a  polynomial  of  degree  2’''ndimn  (ndinm  =  2  or  3) 

To  find  the  point  x  minimizing  <rin  the  multidimensional  space  V  a  Conjugate  Gra¬ 
dient  Algorithm  is  used.  This  algorithm  reduces  the  multidimensional  minimization 
problem  to  a  succession  of  one-dimensional  minimization  problems.  The  different 
steps  are  sketched  below  : 

Initialization  : 

Xo  6  V;  coordinates  of  the  initizii  mesh, 
ho  -  -V<r(xo) 

Assume  that  Xn-i  and  h  are  known, 

1  •  Computation  of  Xn  : 

X„  =  Xn_l  h  n-l 

with,  Pn  =  Arg  mmp<7(x„_i  +  p  h  „_i) 

2  -  Computation  of  7„  : 

with,  G«  =  -V<t(x„) 

3  -  Compute  the  descent  direction  h  : 

bn  =  -i-7n  h 

4  -  convergence  test  : 

if  I  h  nl  ><>”  =  »'+  1 
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Abstract — The  combination  of  adaptive  remeshing  techniques,  flow  solvers  for  transient  problems  with 
moving  grids,  and  consistent  rigid  bixly  motion  integrators  in  three  dimensions  is  presented.  The  resulting 
scheme  allows  the  economical  simulation  of  fully  coupled  fluid-rigid  body  interaction  problems  of 
•s  arbitrary  geometric  complexity.  Several  results,  showing  three-dimensional  store  separation,  are  given  to 

demonstrate  the  capabilities  developed. 


1.  INTRODUCTION 

Practical  military  applications  where  we  encounter  a 
significant  flowfield  perturbation  due  to  body  motion 
are  as  follows;  ordnance  store  separation  from 
aircraft,  interstage  separation  in  rockets,  shroud 
removal  for  interceptors,  separation  of  MIRVs,  and 
torpedo  launch.  Civilian  applications  where  we 
encounter  a  significant  perturbation  of  the  flowfields 
due  to  body  motion  are  as  follows;  reciprocating 
engines,  turbines,  propellers,  ventilators  and  valves. 

1.1.  Store  separation  as  a  model  problem 

In  order  to  focus  the  discussion  and  the  relevant 
ideas,  consider  the  separation  of  ordnance  stores 
from  airplanes  flying  at  supersonic  speeds.  Figure  I 
illustrates  some  of  the  relevant  physical  processes 
involved  in  these  situations. 

— Shock /shock  interactions:  at  supersonic  speeds, 
the  presence  of  shocks  in  the  flowfields  becomes 
unavoidable.  With  several  bodies  interfering  with 
each  other,  the  shocks  emanating  from  them 
interact  with  each  other  in  sometimes  extremely 
complicated  ways.'-^ 

— Shock -boundary  layer  interaction  :  when  shocks 
impact  on  a  surface,  the  boundary  layer  is  greatly 
influenced  by  parameters  such  as  shock-reflection 
angle,  shock-strength,  the  pressure  gradients  up¬ 
stream  and  downstream  of  the  impact  zone,  and 
body  curvature.  The  resulting  flowfield  may 
vary  abruptly  with  only  minor  changes  of  flight 
conditions. 

— Turbulent,  separated  flows:  given  the  highly 
complex,  and  not  aerodynamically  streamlined 
geometries  of  bomb-bays,  many  of  the  flowfields 
contemplated  will  have  vast  regions  of  separated, 
turbulent  flow.  This  implies  that  any  hope  of 


simulating  them  accurately  with  a  conventional, 
algebraic  turbulence  model  has  to  be  forfeited. 
The  lowest  order  turbulence  model  that  can  yield 
acceptable  results  for  flows  of  this  type  is  the  k,  e 
model. 

— Body  motion:  a  futher  degree  of  complexity  is 
added  for  the  class  of  problems  considered 
here  due  to  the  relation  motion  of  the  bodies 
present.  The  bodies  move  through  an  already 
complex,  highly  nonlinear  flowfield,  modifying  it 
constantly. 

All  of  these  aspects,  taken  together,  make  the 
intuitive  prediction  of  these  flows,  as  well  as  any 
extrapolation  from  past  experience,  a  very  unreliable 
design  approach.  The  nonlinear  character  of  these 
flows  also  implies  that  safe  deployment  from  a 
mid-cavity  position  does  not  guarantee  safe  deploy¬ 
ment  from  a  side-cavity  position.  The  only  other  two 
alternative  design  procedures  besides  computational 
fluid  dynamics  (CFD),  wind-tunnel  measurements 
and  flight  testing,  are  either  extremely  expensive  or 
impossible. 

Wind-tunnel  experiments  for  store  separation  in 
the  supersonic  regime  are  difficult  because; 

— Three  non-dimensional  numbers  need  to  be 
reproduced  on  the  scaled  model  at  the  same  time; 
Reynolds  number,  Mach  number  and  Froude 
number.  For  supersonic  flows  in  particular,  the 
reproduction  of  all  three  non-dimensional  numbers 
in  the  wind-tunnel  is  practically  impossible. 

— The  release  of  ordnance  stores  requires  several 
seconds,  a  time-frame  that  would  be  too  power¬ 
consuming — and  thus  expensive — for  most  large 
wind-tunnels. 

— The  release  of  ordnance  stores  into  a  supersonic 
free  stream  tends  to  accelerate  these  objects 
drastically,  propelling  them  to  high  velocities  very 
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Enlargement  of  Area  A 


Fig.  1.  Store  separation:  relevant  physical  processes. 


rapidly.  Thus,  one  can  expect  extensive  damage 
from  any  experimental  program  of  this  sort. 

Therefore,  in-flight  experiments  appear  as  the  only 
viable  choice.  However,  this  is  a  very  primitive,  and 
extremely  expensive,  design  process. 

— A  protoype  has  to  be  built  to  attain  certainty  in 
the  safety  of  the  design.  Production  of  a  complete 
prototype  on  such  uncertain  terms — there  is  no 
guarantee  that  it  will  deploy  safely — appears  almost 
unjustifiable. 

— Unsafe  deployment  may  damage  or  destroy  the 
carrier  vehicle.  This  high  risk  implies  additional 
expenses  in  any  test  program. 

— The  prototype  has  to  be  tested  for  each  new  release 
position;  safe  deployment  from  a  given  position  at 
a  certain  speed  and  height  does  not  imply  safe 


deployment  from  any  other  {.vjsition  and/or  any 
other  speed  and/or  any  other  height.  As  one  can  see, 
this  can  lead  to  an  extremely  lengthy,  and  costly, 
certification  procedure  for  each  new  ordnance  store 
entering  service. 

The  situation  outlined  above  for  store  separation  is 
not  much  diflerent  from  that  encountered  in  all  of 
the  other  engineering  applications  listed  above.  The 
main  difficulty  in  predicting  all  of  these  flowfields  stems 
from  the  fact  that  body  motion  will,  in  most  cases, 
lead  to  complex,  time  dependent  flows.  Advances  in 
computer  speed  and  memory  over  the  next  decade 
will  allow  the  simulation  of  these  flows  on  a  routine 
basis.  Thus,  one  can  expect  CFD  to  gradually  take  the 
lead  role  in  the  design  process  for  these  applications. 
The  present  effort  represents  a  first  step  in  this  direction. 
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1.2.  General  features  of  any  computational  fluid 

dynamics  methodology  for  moving  bodies 

The  accurate  simulation  of  three-dimensional  time- 

dependent,  compressible  flows  with  moving  bodies 

requires  the  following  capabilities. 

(a)  Interactive  grid  generation  methods.  The  fast, 
userfriendly,  interactive  generation  of  grids  in 
three  dimensions  is  essential  to  the  success  and 
widespread  acceptance  of  any  CFD  tool  in  the 
user  community.  Without  such  a  capability,  the 
set-up  times  for  new  problems  run  in  the  order  of 
months,  significantly  reducing  the  benefits  which 
may  be  realized  from  any  CFD  methodology. 
Thus,  interactive  grid  generation  methods  for 
unstructured  grids  that  reduce  set-up  times  to 
days  if  not  hours  are  a  prime  requirement. 

(b)  Solvers  that  can  handle  moving  frames  of  refer¬ 
ence.  Since  at  the  very  least  the  portions  of  the 
mesh  close  to  the  moving  bodies  will  move  in 
time,  the  ability  to  describe  the  equations  of 
motion  for  the  fluids  in  moving  frames  of  reference 
becomes  mandatory. 

(c)  High -order,  monotonicity -preserving  schemes. 
These  schemes  are  needed  to  simulate  time- 
dependent  flows  with  strong  shocks  and  other 
discontinuities  that  will  arise  in  the  flows  of 
interest  in  this  effort  (supersonic  and  hypersonic 
speeds). 

(d)  Proper  modeling  of  turbulent,  separated  flows. 
The  flowfields  considered  will  have  vast  regions 
of  separated,  turbulent  flow.  The  lowest  order 
turbulence  model  that  can  yield  acceptable  results 
for  the  flows  considered  here  is  the  k,  e  model. 

(e)  Fast  regridding  capability  for  regions  with  distorted 
elements.  These  are  needed  because  the  motion 
of  bodies  may  be  severe,  leading  to  distorted 
elements  which  in  turn  lead  to  poor  numerical 
results. 

(f)  Adaptive  refinement  schemes.  Experience  over 
recent  years*^^  has  demonstrated  that  self-adaptive 
refinement  schemes  are  essential  in  reducing  the 
total  number  of  degrees  of  freedom  without 
deteriorating  the  accuracy  of  the  solution.  In  three 
dimensions,  the  ability  to  refine  locally  regions  of 
interest  will  determine  the  accuracy  of  the  result 
and  whether  it  can  be  obtained  in  a  reasonable 
time.  Given  the  currently  available  hardware,  it  is 
impossible  to  solve  three-dimensional  problems 
using  uniformly  fine  grids  everywhere  in  the 
computational  domain. 

(g)  Consistent  rigid  body  motion  integrators.  In  order 
to  fully  couple  the  motion  of  rigid  bodies  with  the 
aerodynamic  forces  exerted  on  them,  consistent 
rigid  body  motion  integrators  must  be  developed. 
This  task  is  relatively  simple  in  two  dimensions. 
However,  in  three  dimensions  the  temporal  vari¬ 


ation  of  the  moments  of  inertia  tensor  can  lead  to 
difficulties. 

(h)  Interactive  post -processing  capabilities.  Under¬ 
standing  of  the  complex,  time-dependent,  three- 
dimensional  flowfields  requires  instantaneous 
visualization  of  several  key  parameters  such  as 
pressure,  Mach  number,  density,  etc.  An  engineer 
that  cannot  visualize  immediately  a  computed 
flowfield,  in  order  to  make  judicious  changes  in 
the  design,  will  never  accept  CFD  as  a  design 
tool.  Thus,  a  fast,  interactive,  workstation-based 
post-processing  capability  is  required. 

This  list  indicates  that  techniques  from  several 
different  areas  of  CFD  and  computer  science  must  be 
combined  to  meet  the  desired  goal.  It  is  therefore  not 
surprising  that  very  few  attempts  have  been  made  to 
tackle  the  complete  class  of  problems.  Currently,  the 
chimera  grid  scheme*  seems  to  be  the  most  promising 
approach  for  structured  grids.  In  this  approach,  local 
grids  for  each  body  are  overset  on  a  major  grid  that 
covers  the  complete  computational  domain.  For  un¬ 
structured  grids,  Formaggia  et  al?  used  local  remesh¬ 
ing  for  regions  of  distorted  elements  in  combination 
with  Eulerian  and  arbitrary  Eulerian-Lagrangian 
solvers  in  two  dimensions  to  simulate  store  separation 
problems.  In  both  cases,  the  body  motion  was  pre¬ 
scribed,  and  no  adaptive  refinement  techniques  were 
employed.  In  1988  the  present  author  presented  a 
fully  coupled  two-dimensional  fluid-rigid  body  inter¬ 
action  algorithm.**  This  algorithm  also  employed 
adaptive  remeshing  to  accurately  simulate  the  flowfield 
at  hand.  This  development  represented  the  first  attempt 
to  combine  and  incorporate  in  a  single,  coherent 
software  package  all  of  the  requirements  listed  above. 

The  present  paper  extends  this  methodology  to 
three  dimensions.  While  conceptually  the  same  as  the 
two-dimensional  algorithm,  the  three-dimensional 
extension  requires  several  important  improvements; 
better  three-dimensional  grid  generators,  consistent 
three-dimensional  rigid  body  motion  integrators, 
interactive  plotting  tools,  and  access  to  a  large  memory 
supercomputer  for  debugging.  Given  the  currently 
available  computer  hardware,  and  our  lack  of  knowl¬ 
edge  in  turbulence  modeling,  it  seems  unreasonable  to 
include  turbulence  modeling  at  the  present  stage  of 
development.  Therefore,  the  present  discussion  will 
center  on  Euler  solvers,  rather  than  Navier-Stokes 
solvers  for  compressible  flows. 

The  rest  of  this  paper  is  divided  as  follows.  Section 
2  treats  the  equations  of  motion  for  the  flowfield  in 
arbitrary  frames  of  reference,  as  well  as  their  solution 
[items  (b)  and  (c)  above].  Section  3  deals  with  the 
equations  of  motion  for  the  moving  bodies  [item  (g)]. 
In  the  present  case,  we  restrict  the  description  to  rigid 
bodies.  Section  4  outlines  the  gridding  technique  used 
[items  (a)  and  (e)].  The  gridding  technique  is  also  used 
to  adaptively  regrid  the  computational  domain  [item 
(f)J.  Finally.  Sec.  5  contains  numerical  examples  that 
demonstrate  the  capabilities  developed. 
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2.  THE  EQUATIONS  OF  MOTION  FOR  THE  FLUID  where  t  an>  ^ne  tangential  and  normal  vectors. 

The  desii  .entum  at  the  new  timestep  should. 

In  order  to  handle  the  moving  frames  of  reference  component  (p  =  0) 

associated  with  the  moving  fimte  elements,  the  partial 
differential  equations  need  to  be  modified.  This  is  most 

easily  accomplished  by  the  Arbitrary  Lagrangian-  Apv"+' =A[p(w  +  at)].  (4) 

Eulerian  (ALE)  formulation.  The  derivation  of  the 

equations  may  be  found  in  Ref  10.  Here,  we  just  state  Combining  Eqs  (3)  and  (4),  we  obtain  for  the  two 

the  final  form  of  the  equations  of  motion.  Given  the  following  cases: 
velocity  field  w  for  the  elements 

w  =  (w\  H-0,  (1)  * 

the  Euler  equations  that  describe  an  inviscid,  com-  Apv"  * '  =  Apw  -t-  [(Apv*  —  Apw)- 1]  •  t;  (5) 

pressible  fluid  may  be  written  as  _ 


p  («'  —  wOp  («''  —  >*'0p 

pu*  («'  —  w^jpu'  +  p  (u’’  —  w’')pu‘ 

pu’’  ►  4-  "  (u*  —  ^  (u”  —  it’’')pa’’  H 

pM*  (u^  —  w")pu^  (u'  —  tcOpa' 

pe  ,  (_(«'  — H'*)pc -f  X 

Observe  that  in  the  case  of  no  element  movement 
(w  =  0),  we  recover  the  usual  Eulerian  conservation- 
law  form  of  the  Euler  equations.  If,  however,  the 
elements  move  with  the  particle  velocity  (w  =  v),  we 
recover  the  Lagrangian  form  of  the  equations  of 
motion.  From  the  numerical  point  of  view,  Eq.  (2) 
implies  that  all  that  is  required  when  going  from  an 
Eulerian  frame  to  an  ALE  frame  is  a  modified 
evaluation  of  the  fluxes  on  the  left-hand  side,  and 
the  additional  evaluation  of  source-terms  on  the 
right-hand  side. 

As  the  elements  move,  their  geometric  parameters 
(shape-function  derivatives,  Jacobians,  etc.)  need  to 
be  recomputed  every  timestep.  If  the  whole  mesh 
is  assumed  to  be  in  motion,  then  these  geometric 
parameters  need  to  be  recomputed  globally.  In  order 
to  save  CPU  time,  only  a  small  number  of  elements 
surrounding  the  bodies  are  actually  moved.  The 
remainder  of  the  field  is  then  treated  in  the  usual 
Eulerian  frame  of  reference,  avoiding  the  need  to 
recompute  geometric  parameters.  This  is  accomplished 
by  identifying  several  layers  of  elements  surrounding 
the  bodies,  which  are  then  moved.  As  the  number  of 
layers  increases,  the  time-interval  between  regridding 
increases,  but  so  also  does  the  cost  per  timestep. 
Therefore,  one  has  to  strike  a  balance  between  the 
CPU  requirements  per  timestep  and  the  CPU  require¬ 
ments  per  regridding.  In  the  present  case,  we  found 
that  two  to  five  layers  of  elements  represented  a  good 
compromise. 

2.1.  Boundary  conditions 

When  imposing  the  boundary  conditions  for  the 
velocities  at  solid  wails,  we  need  to  take  the  velocity 
of  the  surface  w  into  consideration.  Denoting  the 
predicted  momentum  at  the  surface  as  Apv*,  we  can 
decompose  it  as  follows; 


(«'  -  wOp  '  p 

(u^  —  w^)pu^  pu^ 

t- p  ^  (u' —  >V^)pM’’  ►=— V  w-  p«^'  ►. 

(u‘  —  w^)pu^+p  pu‘ 

“'pX  L(“'-‘‘’')pe +u'pJ.,.  Lp^  J 

_ ^ 

(b)  given  n 

Apv"  '  =  Apv*  —  [(Apv*  —  Apw)  •  n]  ■  n.  (6) 

2.2.  The  flow  solver  (FEM-FCT) 

For  the  compressible  flows  described  by  Eqn  (2), 
discontinuities  in  the  variables  may  arise  (e.g.  shocks 
or  contact  discontinuities).  Any  numerical  scheme  of 
order  higher  than  one  will  produce  overshoots  or 
ripples  at  such  discontinuities  (the  so-called  “Godunov 
theorem”).  In  the  present  case  the  appearance  of  these 
overshoots,  which  may  lead  to  numerical  instability, 
is  avoided  by  combining,  in  a  conservative  manner, 
a  high-order  scheme  with  a  low-order  scheme."  The 
temporal  discretization  of  Eq.  (2)  yields 

t/"+'  =  C/"-t-At/,  (7) 

where  AU  is  the  increment  of  the  unknowns  obtained 
for  a  given  scheme  at  time  t  =  t".  Our  aim  is  to  obtain 
a  A  (/  of  as  high  an  order  as  possible  without  intro¬ 
ducing  overshoots.  To  this  end,  we  rewrite  Eq.  (7)  as 

1/"+' =  {/-  + AI/'-KAt/‘-At/')  (8) 

or 

1/-+' =  {/'-(-(At/*- AI/‘).  (9) 

Here  AU*  and  A(/'  denote  the  increments  obtained 
by  some  high-  and  low-order  scheme,  respectively, 
whereas  U'  is  the  monotone,  ripple-free  solution 
at  time  f  =  f "  '  of  the  low-order  scheme.  The  idea 
behind  FCT  is  to  limit  the  second  term  on  the 
right-hand  side  of  Eq.  (9) 


Apv*  =  A[p(w  -t-  art  -(-  pB)), 


(3) 


{/"+'  =  I/'  +  lim(Af/*-At/'), 


(10) 
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in  such  a  way  that  no  new  overshoots  or  undershoots 
are  created.  It  is  at  this  point  that  a  further  constraint, 
given  by  the  conservation  law  (2)  itself,  must  be  taken 
into  account:  strict  conservation  on  the  discrete  level 
should  be  maintained.  The  simplest  way  to  guarantee 
this  for  the  node-centered  schemes  considered  here 
is  by  constructing  schemes  for  which  the  sum  of  the 
contributions  of  each  individual  element  (cell)  to  its 
surrounding  nodes  vanishes.  This  means  that  the 
limiting  process  [Eq.  (10)]  will  have  to  be  carried  out 
in  the  elements  (cells).  Further  deiaiis  uu  the  limiting 
procedure,  its  algorithmic  implementation,  and  the 
high-  and  low-order  schemes  employed  may  be  found 
in  Ref.  11. 

3.  THE  EQUATIONS  OF  MOTION 
FOR  THE  RIGID  BODIES 

The  movement  of  rigid  bodies  can  be  found  in 
standard  textbooks  on  classical  mechanics  (see  e.g. 
Ref.  12).  Due  to  its  nonlinear  character,  rigid  body 
motion  in  three  dimensions  is  not  as  straightforward 
as  it  may  seem.  Therefore,  a  more  detailed  description 
of  the  numerical  implementation  used  is  given  here. 
The  situation  under  consideration  is  shown  in  Fig.  2. 
Given  the  position  vector  of  any  point  of  the  body 

r  =  r,-fro  (H) 

the  velocity  and  acceleration  of  this  point  will  be 
f  =  f,-l-»o  =  »,  +  <oxro  (12) 

?  =  X  ro-t-to  X  (o)  X  fg).  (13) 

Using  the  vector  relationships 

r  X  (to  X  ((B  X  r))  =  (r  •  (oo  X  r))®  —  (r  •  <bX®  x  •■) 

=  —<0  X  (t0r)  ■  0)  (14) 

and  the  abbreviations 


y 


4  =  J  '■o'’ipdi2  (16) 

©  =  fr(I)  •  1  -  I 

^yy  "1”  ^11  ^xy 

-Ixy,  /„  +  /„  -lyzV  (12) 

^xx  ^yz  fjrjc  "1”  ^yy\ 

we  then  have  the  following  equations  describing 
balance  of  forces  and  moments: 

mv,.  =  =  —  J  pndr  (18) 

0©  —  0)  X  (I  •  to)  =  5]ro  X  F  =  —  j"  prgxndr.  (19) 

Observe  that,  in  two  dimensions,  the  second  term  on 
the  left-hand  side  disappears,  considerably  simplify¬ 
ing  the  equations.  However,  in  three  dimensions  it 
usually  does  not.  Another  complication  that  arises 
only  in  three  dimensions  is  the  temporal  variation  of 
the  inertial  matrix  @.  As  one  can  see  from  Eq.  (16), 
the  values  of  ©  will  vary  as  the  body  rotates.  This 
implies  that  during  the  simulation  one  has  to  follow 
the  local  frame  of  reference  of  the  body. 

In  order  to  update  the  velocities  and  positions 
of  the  bodies  in  time,  we  employ  an  explicit  time- 
marching  scheme.  This  seems  reasonable,  as  in 
practical  calculations  the  timescales  of  the  body- 
movement  are  much  larger  than  those  associated  with 
the  fluid  flow.  Thus,  we  update  v^,  ©  as  follows: 


Y^+l  yj-f. 

(20) 

©"+'  =  ©"-i-Ar©". 

(21) 

A  minor  difficulty  now  becomes  apparent:  the  magni¬ 
tude  of  the  timestep  At  is  unknown  before  the  start 
of  the  flowfield  update.  In  the  present  case,  the 
timestep  of  the  previous  timestep  was  taken  instead. 
This  implies  that  the  body  movement  is  “lagging” 
behind  the  flowfield  by  at  most  one  timestep.  How¬ 
ever,  practical  simulations  show  that  the  actual  error 
is  much  smaller,  as  the  magnitude  of  At  does  not 
change  abruptly.  For  the  time-interval  [t",  we 
then  have  the  average  velocities 

C  =  0.5*(v;  +  '-(-v;)  (22) 

©*'’  =  0.5*  (©"+'-!-©”).  (23) 

Combining  Eqs  (22)  and  (23)  with  Eq.  (12),  we  are 
now  in  a  position  to  compute  the  velocities  at  the 
surface  of  the  bodies,  Wp. 

Some  of  the  simulations  shown  below  required 
several  thousand  timesteps.  If  one  simply  uses  the 
velocities  obtained  at  the  boundary  from  Eqs  (22) 
and  (23),  the  body  shape  becomes  more  and  more 
distorted.  This  is  a  purely  numerical  artifact.  It  can 


Fig.  2.  Rigid  body  motion. 
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Fig.  3.  Elongation  of  body.  Correct  path:  AB.  Computed; 
AC. 

be  explained  by  looking  at  the  situation  depicted  in 
Fig.  3.  The  portions  of  the  body  with  higher  velocity 
tend  to  “elongate”  the  body.  This  implies  that  one 
ought  to  impose  the  exact  rigid  body  motion  when 
updating  points  on  the  surface.  With  reference  to 
Fig.  4,  we  decompose  a  point  lying  on  the  body  at 
time  t  —  t"  into  three  components 

r"  =  r,-f  r,  +  r,.  (24) 

We  can  then  define  unit  vectors  in  the  directions  of 
r,  and  r. 


Furthermore,  we  define  the  vector  e,  as 

e,  =  e,  (26) 

Then,  given  the  incremental  rotation  angle  A(p  = 

I  to*  I  At,  the  new  position  for  r  is  obtained  from 

'  =  r,  +  Afv"  +  r,  +  |r,|(cos(A<p)e,  +  sin(A<p)eJ. 

(27) 

The  complete  rigid  body  algorithm  then  consists  of  the 
following  steps. 


Fig.  4.  Decomposition  of  surface  vector  for  rigid  body 
motion. 


(B.l)  Compute  body  forces  and  moments  from 
Eqs  (18)  and  (19). 

(B.2)  Transform  moments  to  the  local  frame  of 
reference  of  the  body 

W'  =  (e'’  •  (28) 

(B.3)  Given  the  estimated  timestep  At,  obtain  the 
accelerations  v^,  <a  from  Eqs  (18)  and  (19). 
(B.4)  Given  the  accelerations,  compute  average 
velocities  vf*,  co**  for  tinie-intei^'al  [f", 
from  Eqs  (22)  and  (23). 

(B.5)  Transform  back  the  angular  velocity  to*  from 
the  local  frame  of  reference  of  the  body  to 
Cartesian  coordinates 

oj'  =  (e'  •  e'')to'\  (29) 

(B.6)  Given  the  actual  timestep  At,  update  the 
positions  of  the  points  lying  on  the  surface  of 
the  body,  as  well  as  the  points  defining  the 
body  geometry  using  Eqs  (24>-(27). 

(B.7)  Given  the  actual  timestep  At,  update  the  pos¬ 
itions  of  the  centers  of  mass  and  the  rotational 
frame  of  reference  using  Eqs  (24)-(27). 

4.  ADAPTIVE  REMESHING 

For  typical  compressible  flow  problems,  we  have 
small  regions  of  rapid  change  in  the  solution  embedded 
in  large  regions  where  the  solution  is  smooth.  In  order 
to  simulate  correctly  the  interaction  of  these  dis¬ 
continuities  or  fronts,  an  appropriately  fine  mesh  is 
required.  It  would,  however,  be  extremely  wasteful 
to  have  an  overall  fine  mesh,  as  the  regions  where  a 
fine  mesh  is  required  are  small.  Therefore,  the  use  of 
adaptive  refinement  techniques  becomes  imperative. 
As  the  bodies  in  the  flowfield  may  undergo  arbitrary 
movement  (see  examples  below),  a  fixed  mesh  structure 
will  lead  to  badly  distorted  elements.  This  means  that 
at  least  a  partial  regeneration  of  the  computational 
domain  is  required.  On  the  other  hand,  as  the  bodies 
move  through  the  flowfield,  the  positions  of  relevant 
flow  features  will  change.  Therefore,  in  most  of  the 
computational  domain  a  new  mesh  distribution  will 
be  required.  The  idea  is  to  regenerate  the  whole 
computational  domain  adaptively,  taking  into  con¬ 
sideration  the  current  flowfield  solution.  In  order  to 
generate  or  regenerate  a  mesh  we  use  the  advancing 
front  technique:^’  ’’ 

(F.l)  Use  the  current  grid  and  solution,  together 
with  appropriate  error  indicators,  to  define  the 
spatial  variation  of  the  size,  the  stretching,  and 
the  stretching  direction  of  the  elements  to  be 
generated.  At  the  nodes  of  the  current  grid  we 
define  the  desired  element  size,  element  stretch¬ 
ing,  and  stretching  direction.  In  what  follows 
we  will  denote  this  grid  as  the  background  grid. 
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(F.2)  Define  the  boundaries  of  the  domain  to  be 
gridded.  This  is  typically  accomplished  by 
splines  in  two  dimensions  and  surface  patches 
in  three  dimensions. 

(F.3)  Using  the  information  stored  on  the  back¬ 
ground  grid,  set  up  faces  on  all  these  bound¬ 
aries.  This  yields  the  initial  front  of  faces.  At 
the  same  time,  find  the  generation  parameters 
(element  size,  element  stretching  and  stretching 
direction)  for  these  faces  from  the  background 
grid 

(F.4)  Select  the  next  face  to  be  deleted  from  the 
front;  in  order  to  avoid  large  elements  crossing 
over  regions  of  small  elements,  the  face  form¬ 
ing  the  smallest  new  element  is  selected  as  the 
next  face  to  be  deleted  from  the  list  of  faces. 
(F.5)  For  the  face  to  be  deleted. 

(F.5.1)  Determine  a  “best  point”  position 
for  the  introduction  of  a  new  point 
IPNEW. 

(F.5. 2)  Determine  whether  a  point  exists  in  the 
already  generated  grid  that  should  be 
used  in  lieu  of  the  new  point.  IF  there 
is  such  a  point,  set  this  point  to  IPNEW 
and  continue  searching  (go  to  F.5. 2). 
(F.5. 3)  Determine  whether  the  element  formed 
with  the  selected  point  IPNEW  does 
not  cross  any  given  faces.  If  it  does, 
select  a  new  point  as  IPNEW  and  try 
again  (go  to  F.5.3). 

(F.6)  Add  the  new  element,  point,  and  faces  to  their 
respective  lists. 

(F.7)  Find  the  generation  parameters  for  the  new 
faces  from  the  background  grid. 

(F.8)  Delete  the  known  faces  from  the  list  of  faces. 
(F.9)  If  there  are  any  faces  left  in  the  front,  go  to  F.4. 

4. 1 .  Recent  developments 

A  typical  simulation  where  bodies  undergo  severe 
motion  typically  requires  several  tens,  if  not  hundreds, 
of  remeshings.  Therefore,  the  grid  generator  must  be 
reliable  and  fast. 

4.1.1.  Reliability.  We  have  recently  increased  the 
reliability  of  the  grid  generator  to  a  point  where  it 
can  be  applied  on  a  routine  basis  in  a  production 
environment.  This  significant  increase  in  reliability 
was  achieved  by:  (a)  not  allowing  any  bad  elements 
during  the  generation  process;  and  (b)  enlarging  and 
remeshing  those  regions  where  new  elements  could 
not  be  introduced.  Thus,  we  first  attempt  to  complete 
the  mesh,  skipping  those  faces  that  do  not  give  rise 
to  good  elements.  If  pockets  of  unmeshed  regions 
remain,  we  enlarge  them  somewhat,  and  regrid  them. 
This  technique  has  proven  extremely  robust  and  reli¬ 
able.  It  has  also  made  smoothing  of  meshes  possible: 
if  elements  with  negative  or  small  Jacobians  appear 
during  smoothing,  these  elements  are  removed.  The 
unmeshed  regions  of  space  are  then  regridded.  By 


being  able  to  smooth,  the  mesh  quality  was  improved 
substantially. 

4.1.2.  Speed.  The  following  means  are  used  to 
achieve  speed. 

(a)  Use  of  optimal  data  structures.  The  operations 
that  could  potentially  reduce  the  efficiency  of  the 
algorithm  to  0(N'  ^)  or  even  0{N^)  are  (see  Sec.  2) 
as  follows. 

— Finding  the  next  fate  to  be  deleted  (step  F.4). 

— Finding  the  closest  given  point  to  a  new  point 
(step  F.5. 2). 

— Finding  the  faces  adjacent  to  a  given  {joint  (step 
F.5.3). 

—Finding  for  any  giving  location  the  values  of 
generation  parameters  from  the  background  grid 
(steps  F.3  ana  F.7).  This  is  an  interpiolation 
problem  on  unstructured  grids. 

The  verb  “find”  apujears  in  all  of  these  o{)erations. 
The  main  task  is  to  design  the  best  data  structures 
for  [jerforming  these  search  o|;)erations  as  efficiently 
as  {wssible.  The  data  structures  used  are  as  follows. 

— Heap-lists  to  find  the  next  face  to  be  deleted  from 
the  from. 

— Quad-trees  (two-dimensional)  and  octrees  (three- 
dimensional)  to  locate  (joints  that  are  close  to  any 
given  location. 

— Linked  lists  to  determine  which  faces  are  adjacent 
to  a  point. 

The  detailed  implementation  of  these  data-structures 
may  be  found  in  Ref.  13. 

(b)  Filtering.  Typically,  the  number  of  close  [joints 
and  faces  is  far  too  conservative,  i.e.  large.  As  an 
example,  consider  the  search  for  close  [joints:  there 
may  be  up  to  eight  points  inside  an  octant,  but  of 
these  only  one  may  be  close  to  the  face  to  be  taken 
out.  The  idea  is  to  filter  out  these  “distant”  faces  and 
points  in  order  to  avoid  extra  work  afterwards.  While 
the  search  operations  are  difficult  to  vectorize,  these 
filtering  operations  lend  themselves  to  vectorization 
in  a  straightforward  way,  leading  to  a  considerable 
overall  reduction  in  CPU  requirements. 

(c)  Automatic  reduction  of  unused  points.  As  the 
front  advances  into  the  domain  and  more  and  more 
tetrahedra  are  generated,  the  number  of  tree-levels 
increases.  This  automatically  implies  an  increase  in 
CPU  time,  as  more  steps  are  required  to  reach  the 
lower  levels  of  the  trees.  In  order  to  reduce  this  CPU 
increase  as  much  as  possible,  all  trees  are  automati¬ 
cally  restructured.  All  points  which  are  completely 
surrounded  by  tetrahedra  are  eliminated  from  the 
trees.  We  have  found  this  procedure  to  be  extremely 
effective.  It  reduces  the  asymptotic  complexity  of 
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the  grid  generator  to  less  than  0{N  log  N).  In  fact, 
in  most  practical  cases  one  observes  a  linear  0(N) 
asymptotic  complexity,  as  CPU  is  traded  between 
subroutine  call  overheads,  and  less  close  faces  on 
average  for  large  problems. 


(30) 


Defining  the  following  “derivative  quantities”: 


(d)  Global  h-refinement.  While  the  basic  advancing 
front  algorithm  is  a  scalar  algorithm,  /(-refinement 
can  be  completely  vectorized.  Therefore,  the  adaptive 
remeshing  process  can  be  made  considerably  faster 
by  first  generating  a  coarser,  but  stretched  mesh, 
and  then  refining  globally  this  first  mesh  w=th  classic 
h-refinement.*  Typical  speed-ups  achieved  by  using 
this  approach  are  I ;  6  to  1:7. 

Currently,  the  advancing  front  algorithm  constructs 
grids  at  a  rate  of  25,000  tetrahedra  per  minute  on  the 
Cray-XMP  or  Cray-2.  With  one  level  of  h-refinement, 
the  rate  is  190,000-200,000  tetrahedra  per  minute. 
This  rate  is  essentially  independent  of  grid-size,  but 
may  decrease  for  very  small  grids. 

4.2.  Local  remeshing 

Practical  simulations  revealed  that  the  appearance 
of  badly  distorted  elements  occurred  at  a  frequency 
that  was  much  higher  than  that  expected  from  the 
element  size  prescribed.  Given  the  relatively  high  cost 
of  global  remeshing,  we  explored  the  idea  of  local 
remeshing  in  the  vicinity  of  the  elements  that  became 
too  distorted.  Thus,  we  proceed  as  follows. 


Z)“  =  , 

^(|(/,.il  +  2-|t/,l-(-|t/,_,|) 

(31) 

D]  =  \ 

+ It/, 

(32) 

D}  =  \ 

(33) 

the  error  on  the  present  (“old”)  grid  is  given  by 


= 


D\  +  D°' 


(34) 


This  implies  that  a  reduction  of  the  current  element 
size  /!'’“  by  a  fraction  ^  to 


/i“"* 


=  4  • 


(35) 


will  lead  to  the  following  estimated  errors: 

'  D]\+Dr 


(36) 


Thus,  given  the  desired  error  value  £”“,  the  reduction 
factor  ^  becomes 


(L.l)  Identify  the  badly  distorted  elements  in  the 
layers  that  move,  writing  them  into  a  list 
LEREM(1:NEREM). 

(L.2)  Add  to  this  list  the  elements  surrounding  these 
badly  distorted  elements. 

(L.3)  Form  “holes”  in  the  present  mesh  by: 

(L.3.1)  Forming  a  new  background  mesh  with 
the  elements  stored  in  the  list  LEREM. 
(L.3. 2)  Deleting  the  elements  stored  in  LEREM 
from  the  current  mesh. 

(L.3. 3)  Removing  all  unused  points  from  the 
grid  thus  obtained. 

(L.4)  Recompute  the  error  indicators  and  new 
element  distribution  for  the  background  grid. 
(L.5)  Regrid  the  “holes”  using  the  advancing  front 
method. 

Typically,  only  a  very  small  number  of  elements 
( <  10)  becomes  so  distorted  that  a  remeshing  is 
required.  Thus,  local  remeshing  is  a  very  economical 
tool  that  has  allowed  us  to  reduce  CPU  requirements 
by  more  than  60%  for  typical  runs. 

4.3.  Determination  of  element  sizes 

In  order  to  estimate  the  element  size,  stretchings, 
and  stretching  directions,  we  employ  the  modified 
interpolation  theory  error  indicator  proposed  in  Ref  3. 
In  one  dimension,  on  a  uniform  grid  of  element  size 
h.  this  error  indicator  reduces  to  the  following  form: 


£Old  2 


/>,'  + 


'((i>,' 


£<M 

+  — [D,!  + 


on) 


W  +  Dl] 


(37) 


Notice  that  if  the  solution  is  smooth,  implying 
D'  4  />“,  then  the  reduction  factor  reverts  to 


consistent  with  the  second-order  accuracy  assumption 
of  linear  elements.  However,  close  to  a  discontinuity, 
where  D'pD'^,  the  reduction  factor  f  is  given  by 

^oew 

(39) 


In  two  and  three  dimensions  we  define  the  corre¬ 
sponding  matrices 


(D%  =  /.V, 

,  [  |.vi|lNi||u,ldn 

(40) 

Jo 

(DX  =  h^ 

|A/i|lNit;,ldn 

(41) 

• 

(t>%  =  6^| 

j’^/vi/v^dnu,|. 

(42) 
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where  N'  denotes  shape-function  of  node  /,  and  h  is 
a  typical  element  length.  Given  these  matrices,  we 
obtain  the  error-indicator  matrix  E  and  its  modal 
decomposition 


£,v.. 

0 

E,y 

0 

E22 

0  [x 

0 

0 

£33] 

Each  principal  direction  is  then  treated  as  a  one¬ 
dimensional  problem.  Using  Eq.  (37),  we  obtain  three 
diilcreni  element  sizes  S,,  o^,  along  the  principal 
directions.  This  information  is  then  used  to  regenerate 
a  better  grid  for  the  problem  at  hand.  We  remark  on 
the  following  characteristics  of  the  present  error 
indicator. 

(a)  The  error  indicator  is  non-dimensional.  Therefore, 
several  variables  may  be  monitored  at  the  same 
time  in  order  to  accurately  track  all  physical 
phenomena  present.  Thus,  we  can  monitor  both 
density  (shocks,  contact  discontinuities)  and 
vorticity  (boundary  layers)  for  viscous  flow 
problems. 

(b)  The  error  indicator  is  bounded.  This  implies  that 
the  user  does  not  have  to  change  specified  error 
tolerances  from  run  to  run.  We  have  found  that 
for  large  classes  of  problems  the  specified  error 
tolerances  could  be  left  untouched  without 
impediment  to  the  adaptation  process.  We  find 
this  of  particular  value  for  the  non-expert  user 
environment. 

Before  proceeding  to  the  overall  algorithm,  we 
summarize  the  steps  required  for  one  adaptive 
remeshing  as  follows. 

(R.l)  Obtain  the  error  indicator  matrix  for  the  grid- 
points  of  the  present  grid. 

(R.2)  Given  the  error  indicator  matrix,  obtain  the 
element  size,  element  stretching  and  stretching 
direction  for  the  new  grid. 

(R.3)  Using  the  old  grid  as  the  “background  grid,” 
remesh  the  computational  domain  using  the 
advancing  front  technique. 

(R.4)  If  further  levels  of  global  A -refinement  are 
desired:  refine  the  new  grid  globally. 

(R.5)  Interpolate  the  solution  from  the  old  grid  to  the 
new  one. 

5.  THE  OVERALL  ALGORITHM 

The  overall  algorithm  for  the  advancement  of  the 
solution  in  time  is  as  follows. 

(0.1)  Advance  the  solution  one  timestep. 

(A.l)  Compute  the  body  forces  and  moments 
from  the  pressure  field  and  any  exterior 
forces. 


(A. 2)  Taking  into  consideration  the  kinematic 
constraints  for  the  body  movements, 
update  the  velocities  of  the  bodies  at 
1=1"*':  0)"^'.  At  the  same  time, 

obtain  the  average  velocities  for 

the  time-interval  [r", 

(A. 3)  With  the  average  velocities  co",  obtain 
the  velocities  on  the  surface  of  each 
body  for  the  time-interval  [n,  t"*']. 

(A. 4)  Given  the  surface  velocities  Wp  on  the 
boundaries  of  the  global  domain,  obtain 
the  global  velocity  field  for  the  element 
movement. 

(A. 5)  Advance  the  solution  by  one  timestep 
using  the  ALE-FEM-FCT  solver.  This 
yields  the  actual  timestep  At", 

(A.6)  Given  the  actual  timestep  At"  and  the 
velocity  field  for  the  element  movement 
Wn,  update  the  coordinates  of  the  points. 
(A. 7)  Update  the  shape-function  derivatives 
and  other  geometric  parameters  for  the 
elements  that  have  been  moved. 

(A.8)  Update  the  centers  of  mass  r,  for  the 
bodies,  as  well  as  the  coordinates  of  the 
points  defining  the  body  geometry. 

(0.2)  If  the  grid  has  become  too  distorted  close  to  the 
moving  bodies;  adaptively  remesh  these  regions. 
(0.3)  If  the  desired  number  of  timesteps  between 
global  remeshings  has  elapsed:  adaptively 
remesh  the  complete  computational  domain. 
(0.4)  If  the  desired  time-interval  has  elapsed:  stop. 
Otherwise,  advance  the  solution  further  (go  to 
0.1), 

6.  NUMERICAL  EXAMPLES 

We  consider  two  numerical  examples  that  demon¬ 
strate  the  effectiveness  of  the  algorithms  developed. 
In  both  cases  an  idealized  store  release  from  a  bay  at 
supersonic  speed  (Ma^  =  2.0)  is  simulated.  Because 
of  symmetry,  only  half  the  flowfield  domain  needs  to 
be  simulated.  Release  into  a  suptersonic  flowfield  will 
necessitate  the  forceful  ejection  of  stores.  Therefore, 
the  motion  of  the  stores  was  prescribed,  and  the 
resulting  forces  computed.  Adaptive  remeshing  was 
performed  every  100  timesteps  initially,  while  at  latter 
times  the  grid  was  modified  every  40  timesteps. 
The  maximum  stretching  ratio  specified  was  S  =  1.5. 
Density  and  the  absolute  value  of  the  velocity  were 
chosen  as  indicator  variables.  The  latter  provided 
suitable  mesh  adaptation  to  the  shear  layers  in  the 
cavity.  Even  though  the  grid  shows  considerable 
variation  in  element  size,  the  average  grid  size  was  of 
the  order  of  280,000  tetrahedra  for  the  first  example, 
and  350,000  tetrahedra  for  the  second  one.  On  a 
uniform  mesh,  the  required  number  of  elements 
would  have  increased  by  more  than  an  order  of 
magnitude.  The  required  CPU  time  for  runs  of  this 
kind  is  of  the  order  of  several  CRAY-XMP  processor 
hours. 
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Fig.  5.3.  Single  store  separation.  Surface  pressure  at  T  =  0.0. 


6.1.  Single  object  falling  into  supersonic  free  stream 

The  computational  domain  is  shown  in  Figs  S.l 
and  S.2.  Observe  that  the  doors  of  the  bay  simulated 
are  somewhat  thicker  than  in  real  life.  The  store  is  a 
slender  object,  resembling  a  missile.  Figures  5.1  and 
5.2  show  the  mesh  on  the  surface  of  the  compu¬ 
tational  domain  at  time  T  =  0.0.  The  corresponding 
pressure  contours  (30)  are  shown  in  Fig.  5.3.  Figures 
5.4-5. 6  show  the  surface  mesh  and  the  pressure 
contours  (60)  at  time  T  =  8.5.  One  can  clearly  discern 
the  extent  of  mesh  adaptation,  as  well  as  the  consider¬ 
able  change  in  shock-strengths  and  shock-positions 
that  occurred  due  to  body  motion. 

6.2.  Multiple  objects  falling  into  supersonic  free  stream 

The  computational  domain  is  the  same  as  before. 
The  two  stores  resemble  bombs.  The  store  at  the  back 
of  the  cavity  is  ejected  first,  followed  by  the  store  in 
the  front  of  the  cavity.  Figures  6.1  and  6.2  show  the 
mesh  on  the  surface  of  the  computational  domain  at 
time  T  =  0.0.  The  corresponding  pressure  contours 
(60)  are  shown  in  Fig.  6.3.  Figures  6.4-6.6  show  the 
surface  mesh  and  the  pressure  contours  (60)  at  time 
T  =  2.62.  While  the  store  at  the  back  of  the  cavity  has 


already  moved  a  considerable  distance,  the  store  in 
the  front  has  just  begun  to  move.  As  before,  one 
can  clearly  discern  the  extent  of  mesh  adaptation,  as 
well  as  the  considerable  change  in  shock-strengths  and 
shock-positions  that  occurred  due  to  body  motion. 

7.  CONCLUSIONS 

We  have  demonstrated  how  the  combination 
of  adaptive  remeshing  techniques,  flow  solvers  for 
transient  problems  with  moving  grids,  and  integrators 
for  rigid  body  motion  allows  the  simulation  of  fully 
coupled  fluid-rigid  body  interaction  problems  of 
arbitrary  geometric  complexity  in  three  dimensions. 
The  overall  reduction  in  CPU  times  as  compared 
to  those  of  uniformly  fine  grids  depends  strongly  on 
the  stretching  ratios  allowed  by  the  physics  of  the 
problem,  but  typically  lies  between  10  and  50.  Areas 
that  deserve  further  study  are  as  follows: 

— the  diffusive  effect  of  interpolation  while  remeshing 

— extension  to  Navier-Stokes  problems 

— treatment  of  multifluid  interactions 

— extension  to  flexible  bodies  and  structures. 


Fig.  6.4.  Multiple  store  separation.  Surface  mesh  at  T  =  2.62. 


Fig.  6.5.  Multiple  store  separation.  Surface  mesh  at  F  =  2.62. 
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Fig.  6.6.  Multiple  store  separation.  Surface  pressure  at  7"  =  2.62. 
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